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INTRODUCTION 


• * 


Accurate  pointing  information  to  selected  stars  is  required  by  various 
SAMTEC  instrumentation  for  data  acquisition  and  instrumentation  calibration. 

Until  now,  the  program  STASHO  on  the  709'’  computer  was  used  for  this  purpose. 

While  valuable  when  originally  designed,  new  algorithms  and  programming  tech- 
niques developed  by  the  Naval  Observatory  have  improved  the  accuracy  obtainable 
by  programs  like  STASHO.  Operationally,  the  large  amount  of  input  information 
required  for  STASHO  runs  left  a wide  margin  for  the  occurrence  of  human  error. 

For  these  reasons  it  was  decided  to  design  an  observed  position  star  look 
angle  program  based  on  the  mathematics  produced  by  the  Naval  Observatory 
(reference  1).  This  work  was  completed  under  the  ROTI  Capability  Task  (A713b) 
and  under  the  Data  Product  Line  Enhancement  Task  {A7102).  It  was  decided  to 
include  precession-nutation  and  Vernal  Equinox  data  internally  within  the 
program.  An  accuracy  design  goal  of  0.01  or  less  arcseconds  for  the  non-refracted 
pointing  information  was  selected.  Accordingly,  all  physical  effects  (except 
refraction)  that  can  contribute  more  than  0.005  arcseconds  are  corrected.  Since 
it  is  preferable  to  remove  refraction  errors  from  measured  pointing  angles  by 
use  of  Rawinsonde  data  before  regression,  it  was  decided  to  use  only  a minimal 
refraction  model  in  this  program.  However,  since  many  models  of  refraction 
effects  are  available  at  WTR,  this  module  was  programmed  so  that  additional 
refraction  techniques  can  easily  be  added  by  the  user  as  subroutines. 

As  work  in  the  STAR  program  progressed,  it  became  evident  that  additional 
programs  were  needed  to  supplement  the  original  observed  position  star  look 
angle  program.  Tiie  complete  list  of  programs  is; 

STARl  - Tnis  program  takes  mean  star  coordinates  referenced  to 
the  nearest  Bessel ian  New  Year  and  computes  azimuth  and 
elevation  data  for  a given  site  during  a given  time  period. 

STAR2  - This  program  takes  mean  star  coordinates  for  a given  epoch 

and  transforms  them  into  mean  star  coordinates  at  a new  epoch. 

STaR3  - This  program  computes  i/5  of  the  yearly  data  deck  for  STARl 
for  any  given  year  and  outputs  the  data  on  punched  cards. 
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JJJPP  - This  program  takes  a new  BLOCK  DATA  (for  use  in  STAR!) 

and  produces  tables  to  check  against  the  tables  available 
i n the  Anierican  Ephemeris  and  Nautical  Almanac  for  the 
year  in  question. 

Although  the  other  programs  have  been  checked  out,  this  report  will  be 
restricted  to  the  discussion  and  verification  of  STARl. 

The  other  programs  were  much  easier  to  verify,  in  any  case.  Moreover,  a listing 
of  JJJPP  can  be  found  in  Appendix  VII.  The  operating  instructions  for  STAR2 
can  be  found  in  Appendix  II.  For  those  who  wish  the  output  of  STARS,  please 
contact  Dr.  Sinclair  at  FEC/Performance  Analysis  Department.  Also,  look  at 
Appendix  V,  which  concerns  use  of  the  output  of  STARS. 

1 .1  Outline  of  Report 

Section  2 is  devoted  to  a discussion  of  the  processing  and  mathematics 
in  the  STARl  MAIN  program.  The  mathematics  and  processing  involved  in  the 
various  subroutines  are  outlined  in  Section  S.  Numerical  tests  and  numerical 
examples  of  the  procedures  unique  to  STARl  as  compared  to  STASHO  (the  preset^* 
star  look  angle  program  used  at  WTR)  are  presented  in  Section  4.  Also  in- 
cluded in  Section  4 is  a comparison  of  the  various  types  of  Polar  Motion 
data.  The  many  subroutine  outputs  are  verified  in  Section  5.  Finally,  the 
complete  program  verification  and  the  report  conclusions  can  be  found  in 
Section  6.  The  technical  details  of  running  the  program  (operating  instruc- 
tions, flow  diagrams,  listings,  sample  output,  etc.)  are  relegated  to  the 
appendices. 
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2.0 


MAIN  PROGRAM  DESCRIPTION  AND  MATHEMATICS 


The  STAR  1 program  consists  of  a main  program  and  nine  subroutines.  The 
routines  and  their  functions  are; 

. MAIN  Handles  all  I/O,  initial  and  time  dependent 

preprocessing,  and  the  primary  mathematics 
required  for  updating  star  positions. 

. ASDCRK  This  subroutine  transforms  the  Right  Ascension, 

declination  coordinate  system  to  cartesian  coordi- 
nates or  vice  versa. 

. CNSTNT  This  subroutine  updates  the  earth's  obliquity,  the 

equation  of  equinoxes  and  its  first  derivative,  and  the 
secular  changes  in  the  earth's  rotation  rate. 

It  also  computes  the  first  and  second  differences 
for  the  updated  Bessel ian  constants  A,  B,  C,  D, 
and  Independent  constant  f. 

CNSTJ  This  subroutine  is  an  additional  entry  within 

CNSTNT,  It  computes  second-order  corrections 
on  a ten-day  basis. 


This  subroutine  converts  degrees  to  site  octal 
output.  It  has  a second  entry,  OCTO,  which 
sets  the  initial  parameters  for  OCT  with  respect 
to  the  encoder  bit  size  of  the  site. 

This  is  a double  precision  subroutine  which 
premultiplies  a 3 x 1 vector  by  a 3 x 3 matrix. 


This  subroutine  changes  integration  time  to  UTC  time 
and  to  local  time  in  days,  hours,  minutes,  and  decimal 
minutes . 
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^ . VERNAL  This  subroutine  computes  the  mean  longitude  of  the 

Vernal  Equinox  at  midnight,  UTC. 

. BLOCK  DATA  This  subroutine  supplies  the  Bessel ian  day  numbers, 
yearly  constants,  geodetic  constants,  and  initial 
values  for  the  program. 

2.1  MAIN  Program 

The  following  is  a description  of  the  pertinent  processing  flow  through  the 
MAIN  program.  Detailed  descriptions  of  counter  loops  and  miscellaneous 
computations  will  be  addressed  in  Appendix  III. 

A.  Initial  Input 

General  input  to  the  program  is: 

. Number  of  stations,  number  of  stars  per  station 
. Start  day,  stop  day 
. Start  hour  and  stop  hour  for  each  day 
. Prediction  increment  size 

. Station  geodetic  and  astronomic  coordinates 

. Special  station  parameters 

4* 

. Polar  motion  parameters 
. Mean  coordinates  of  the  stars 

First,  the  station  coordinates  are  read  and  converted  to  radians  in  a 
right-handed  system.  Then  the  site-dependent  part  of  the  latitude  rotation 
is  computed: 

QM  (2,2)  = DSIN  (ALAT)  = QM  (3,3) 

QM  (2,3)  = DCOS  (ALAT)  = -QM  (3,2) 

Where  ALAT  is  the  astronomic  latitude  of  the  site. 

I Polar  motion  corrections  are  usually  omitted  in  prediction  runs. 
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Then  the  diurnal  aberration  correction  factor  is  computed. 


where 


AEL  = DSQRT(1.  - E^  DSIN^(DLAT) ) , and 
DIURN  = 0.3198  * (AEL  + HT)  * DCOS(DLAT)/ FACTOR, 


AEL  = radius  of  earth  at  the  site  location  (measured  in  units  of 
Earth  semimajor  axis), 

OIURN  = diurnal  correction  constant, 

E = eccentricity  of  Earth, 

HT  = Height  of  site  above  ellipsoid  (measured  in  units  of  Earth 
semimajor  axis), 

DLAT  - geodetic  latitude  of  site,  and 

FACTOR  = conversion  factor  to  convert  arcseconds  to  radians. 


Finally,  in  the  call  to  OCTO,  the  subroutine  OCT  is  initialized  by  the  bit 
size  of  the  site  encoder. 

If  polar  motion  corrections  are  used,  they  are  read  in  at  this  time.  The 
input  is: 

DUTDOT  = First  derivative  of  P0L(1,  ), 

P0L(1,  ) = Polar  motion  time  correction  (one  value  for  each  day 

of  run) , 

P0L(2,  ) = X-value  of  polar  motion  (one  value  for  each  day  of 

run),  and 

P0L(3,  ) = Y-value  of  polar  motion  (one  value  for  each  day  of  run). 

The  initial  run  parameters  (start  day,  stop  day,  start  hour,  stop  hour),  are 
converted  to  UTC  time  here.  If  any  of  the  following  conditions  hold,  the  run 
is  terminated  with  an  appropriate  printout: 

a.  If  start  day  is  a zero, 

b.  If  stop  day  is  later  than  January  2 of  the  next  year,  and 


c.  If  July  1 falls  betv^en  start  day  and  stop  day. 


The  reason  for  tests  1 and  2 is  that  the  Bessel ian  day  constants  only  cover  one 
year.  The  reason  behind  test  3 is  that  the  Bessel ian  day  constants  are  referenced 
to  the  beginning  of  the  nearest  Bessel ian  year  (about  January  1).  Hence,  they 
are  discontinuous  on  or  about  July  1. 

B.  Daily  Corrections  and  Star  Coordinates 

At  this  point  the  star  coordinates  are  read  in.  They  are: 

. mean  Right  Ascension,  epoch  at  nearest  beginning  of  a 
Bessel ian  year, 

. mean  declination,  same  epoch  as  Right  Ascension, 
absolute  parallax, 

. proper  motion  in  declination,  and 
. proper  motion  in  Right  Ascension. 

The  proper  motion  parameters  may  be  set  to  zero  with  little  loss  of  accuracy. 

The  daily  parameters  are  computed  now.  If  the  polar  motion  corrections  are 
used,  their  effect  is  also  computed  here: 

DAZ  = (P0L(2,  )*DSIN(DLONG)  + P0L(3,  )*DC0S(DL0NG) )/DC0S(DLAT) 

TIM  = P0L(1 , ) 


where 


DAZ  = azimuth  correction, 

P0L(1,  ),  P0L(2,  ),  P0L(3,  ) and  DLAT  have  been  previously  defined, 
DLONG  = geodetic  longitude  of  the  site,  and 
TIM  = time  correction  from  UTC  time  to  UTl  time. 

Provision  was  made  at  this  point  to  read  in  a daily  atmosphere  for  use  in 
sophisticated  refraction  subroutines. 
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The  next  corrections  or.cur  in  a call  to  the  subroutine  (.N‘>lNi.  Ihe 
of  this  subroutine  are  the  Earth's  oblicjuity  (dihedral  angle  between  Lh(!  plarn; 
of  the  equator  and  the  plane  of  the  Earth's  orbit  around  the  sun),  and  the  annual 
precession  in  declination  - both  referenced  to  the  begintiing  of  nearest  hesselian 
year  - and  the  equation  of  equinoxes  and  its  first  derivative,  the  secular  change 
in  the  rotation  rate  of  the  Earth,  and  the  Besselian  day  numbers  A,  15,  C,  and  D 
along  with  the  Independent  day  number  f,  together  with  their  first  and  second  diff- 
erences, all  referet\ced  to  0.0  hour  ephemeris  time  (ET)  of  the  present  day.  The 
equations  for  these  parameters  will  be  discussed  in  the  section  on  the  subroutine 
CNSTNT. 

Initially,  the  Right  Ascension  and  declination  are  corrected  for  propei  motion 
and  transformed  to  the  X,  Y,  Z direction  cosines  (called  XV(1),  X\/(?),  and  XV(3), 
respectively).  This  is  a right-handed  coordinate  system  with  the  X-Y  plane 
coincident  with  the  mean  equator  and  the  positive  X-axis  directed  along  the  mean 
equinox  of  the  nearest  beginning  of  the  Besselian  year.  The  next  cc'rrci  timi,  h 
second  order  correction,  is  computed  every  ten  days,  starting  with  the  firsi,  day 
of  the  run.  These  corrections  are: 

DECO  = DEC  + UD  * TTAU 

RASO  = RAS  + UR  * TTAU 

and 

D2DEC  - XJP  * DTAN  (ULtO) 

D2RAS  = XJJ  * DTAN  (DECO)  * DTAN  (DECO) 

where  DEC  is  the  mean  deilination,  UO  is  the  proper  motion  in  declina.ion.  RAS 
is  the  mean  Right  Ascension,  UR  is  the  proper  motion  in  Right  Ascension.  IIAU 
is  the  fraction  of  the  tropical  year  which  has  elapsed  from  the  begintjing  of 
the  nearest  Besselian  year,  D2DEC  is  the  second  order  correction  in  declination, 
D2RAS  is  the  second  order  correction  in  Right  Ascension,  and  XJP  and  XJJ  are 
constants  obtained  from  subroutine  CNSTJ,  a second  entry  into  CNSTNT.  The 
formulas  for  XJP  and  XJJ  will  be  displayed  in  the  section  for  the  subvoutine 
CNSTNT. 
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The  final  parameters  determined  at  this  time  are  the  East  longitude  of  the 
Vernal  Equinox  at  0.0  hour  UTl  and  the  Earth  spin  rate  with  respect  to  the 
moving  Vernal  Equinox.  The  rate  is: 

WED  = (WE  - DW)  (1.  + DUTDOT).  where 
WED  is  the  current  spin  rate, 

WE  is  the  spin  rate  at  epoch  1900.0, 

DW  is  the  secular  decrease  in  the  spin  rate  as  computed  in  CNSTNT,  and 
DUTDOT  has  been  previously  defined. 

C.  Diurnal  Parallax  Corrections 

For  those  stars  with  a parallax  greater  than  0.01  in  absolute  value,  the  following 
corrections  to  the  X,  Y,  Z coordinates  are  performed. 

XI  = -AA(3)  * PAR/(CABER  * DCOS(E))  - X axis 
X2  = -AA(4)  * PAR  * DCOS(E)/CABER  - Y axis 
X3  = X2  * DTAN(E)  - Z axis 

where 

PAR  is  the  parallax  for  this  star, 

F is  the  Earth's  instantaneous  true  obliquity  (as  of  0.0  hour,  ET), 

CABER  = 20.496  arcseconds,  and 

AA(3)  and  AA(4)  are  the  Besselian  constants  C and  D,  respectively, 
updated  to  0.0  hour,  ET,  of  the  present  day. 

D.  Mean-to-Apparent  Transformation 

The  Besselian  constants  are  updated  by  use  of  first  and  second  differences  via 
the  following  formulae. 

XLl  = (LT  + EPH)/8.64D4 
IF(XL1.GT.0.5)  HLl  = 1.  - XLl 


i 

I 

I 
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where 


C(I)  ^ AA(l)  + OA(I)  * XLl  + 02A{I)  * XLl  * XLl/2.  I - 1.5 
IF  (XLl. GT. 0.5), 

C(I)  = AAP(I)  - DAP(I)  * HLI  + D2AP(I)  * HL1  * HLl/2.  I = 1,5 


LT  = time  passed  since  0.0  hour,  UTC,  in  seconds, 

EPH  = ET  - UTC  in  seconds, 

C(I)  = updated  Bessel ian  constants, 

AA(I),  Besselian  constants  at  0.0  hour,  ET, 

AAP(I)  = same  constants,  next  day, 

DA(I)  and  D2A(I)  = first  and  second  differences  for  this  day,  and 
DAP(I)  and  D2AP(I)  = first  and  second  differences  for  the  next  day. 

Note  that  the  Besselian  constants  are  interpolated  by  first  and  second  differences 
for  + 1/2  day  about  0.0  hour  ET. 

The  next  effect  is  that  of  annual  aberration: 

XVl(l)  = XV(1)  - C(4) 

XV1(2)  = XV(2)  + C(3) 

XVI (3)  = XV(3)  + C(3)  * DTAN(E) 

where  C(3)  aitd  C(4)  are  tne  updated  Besselian  constants  C and  D,  respectively. 

The  other  parameters  have  been  previously  identified.  The  precess lon-nutation 
effects  are  accounted  for  by: 


■'  r T 

1 '^XV2(1)| 

1. 

-C(5) 

[-C(l)  - C(2)  * C(5)]‘ 

’ XVl(l) 

j ; V2 

j j 1 

C(5) 

1. 

[C(2)  - C(l)  * Cvj)] 

XVI (2) 

j [ XV2(3)J 

C(l) 

-C(2) 

1. 

XV1(3) 
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The  next  effect  is  the  addition  of  the  parallax  effect  previously  calculated; 

XVI (1)  = XV2(1)  + XI 
XVI (2)  = XV2(2;  + X2 
XVI  (3)  --  XV2(3)  + X3 

The  coordinates  are  transformed  to  the  Right  Ascension  - declination  system 
and  corrected  for  second-order  effects  previously  calculated. 

Call  ASDCRX  (XVl,  PASl , DECl.l) 

RAS1  -=  RASl  r D2RAS 
DECl  = OFCI  + 02DLC 


The  pair  RASl , DECl  now  represent  the  apparent  position  of  the  star, 

E.  Apparent-to-Observed  Transformations 

First,  the  apparent  Vernal  Equinox  and  hour  angle  are  computed. 


XLll  = LT/8.64D4 
E0FE2  = E0!"E  + XLll  * DEFE 
EI.NG  = ELNGO  + E0FE2 
TIMl  = TIM 

TIMl  = TIMl  + OTDOT  * 10-3  * XLll  if  TIM  f 0. 
OUT  = WED  * (LT  + TIMl) 

ANGL  = ELNG  - DUT 

HA  = 2.  * PI  - (RASl  + ANGL) 


iO 


where 


LT  is  the  time  from  0.0  hour,  UTC, 


EOFE  is  the  equation  of  equinoxes  (difference  between  apparent 
and  mean  equinox)  in  radians. 


DEFE  is  the  change  between  EOFE  now  and  EOFE  tomorrow, 

TIM  is  the  difference  between  UTl  and  UTC  at  0.0  hour,  UTC, 

OTDOT  is  the  rate  of  change  of  UTl  - UTC  in  milliseconds, 

TIMl  is  the  updated  TIM, 

E0FE2  is  the  updated  equation  of  equinoxes, 

ELNG  is  the  updated  apparent  longitude  of  the  Vernal  Equinox  without 


Earth  spin, 

OUT  is  the  change  in  the  longitude  of  the  Vernal  Equinox  caused  by 


Earth  spin  since  0.0  hour,  UTl, 


ANGL  is  the  updated  apparent  longitude  of  the  Vernal  Equinox, 
referenced  to  updated  UTl  (if  TIM  ^0.),  and 


HA  is  the  hour  angle  of  the  star. 


ihe  diurnal  aberration  effect  can  nc.v  be  reiixived: 


RASl  = RASl  + DIURN  * DC0S(HA)/DC0S(DEC1 ) 
DEC!  = OrCl  -r  DIURN  * DSIN(HA)  * DSIN(DECl) 


/here  all  symbols  have  .-.'iready  been  .aentified 


RASl  anc  CltCl  are  then  transformed  to  X,  Y,  Z coordinates. 


r 


The  instantaneous  Right  Ascensiott  of  t^e  site  is  then:  RA  = ALONG  - ANGL, 
where  ALONG  is  the  astronomic  lonoHui^o,  thus  correcting  for  deflection  of 
vertical.  The  star  oosition  is  then  rotated  to  the  site  position  in  longi- 
tude and  latitude  by: 


XV(1) 

1.  0.  0.  ' 

'-DSIN(RA)  DCOS(RA)  0. 

XVl(l) 

XV(2) 

= 

0.  DSIN(ALA^)  DCOS(ALAT) 

-DC0S{RA)  -DSIN(RA)  0. 

XV1(2) 

XV(3) 

0.  -DCOS(DLAT)  DSIN(AIAT) 

0.  1., 

_ XV1(3) 

Above,  ALAT  is  the  site  astronomic  latitude.  Since  geocentric  aberration  can 
be  ignored  for  stars,  there  is  no  trans'formation  along  the  Earth's  radius  to 
the  site.  Hence: 

DLL  = DSQRT(XV(1)'*XV(1)  + XV(2)*XV(2)) 

EL  = DATAN(XV(3)/DEL) 

One  then  checks  to  see  that  the  EL  is  larger  than  the  minimum  EL  for  this  site. 
Otherwise,  the  process  restarts  at  LT  INC  seconds  later. 

The  next  calculation  is  AZ: 

AZ  = DATAN2(XV1(1),  XVI (2))  * RAD  + DAZ 

Above,  AZ  is  azimuth,  EL  is  elevation,  and  DAZ  is  the  polar  motion  correction 
to  azimuth. 

Now,  if  desired,  the  elevation  can  be  corrected  for  refraction.  The  simple 
correction  now  used  is  listed  below.  This  model  is  surprisingly  accurate 
above  12°  elevation,  but  loses  accuracy  at  low  elevation  angles.  However,  a 
more  sophisticated  refraction  subroutine(s)  may  easily  be  added  at  this  point. 

EL  = EL  + 3.36D-A  * RAD/DTAN ( EL/RAD ) . 

The  computationally  faster  "EL  - DARSIN(XV(3))"  was  found  to  be 
inadequate  due  to  the  inaccuracy  of  the  DARSIN  subroutine  when  compared 
to  the  DATAN  subroutine.  See  Section  5 for  a discussion  of  this  difference. 
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Plunge  and  mil  outputs  are  computed  via  the  following  equations. 


I 

i 


I 


I 
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AZMIL  = AZ  * DMIL 

ELMIL  = EL  * DMIL 

ELPLG  = (180.  - EL)  * DMIL 

AZPLG  = DMOD  ((AZ  + 180.),  360.)  * DMIL 

The  octal  equivalents  of  AZ,  EL,  AZPLG,  and  ELPLG  are  computed  via  a call  to 
OCT. 

At  this  time  the  printout  times  are  updated  via: 

XMG  = XMG  + DINC/60D0 
XML  = XML  + DINC/60D0 

IF((XMG.GE.60.).0R.(XML.GE.60.).0R.(NQ.EQ.1))  CALL  TIMEE 

where 

XMG  = number  of  minutes  after  the  UTC  hour, 

XML  = number  of  minutes  after  the  local  time  zone  hour, 

OINC  = integration  step  size,  and 

TIMEE  is  the  subroutine  which  converts  the  variable  LT  to  days,  hours, 
minutes  ..nd  decimal  minutes  in  GMT  and  local  time  for  printout  pur- 
poses. YiMEE  will  be  discussed  in  Section  3. 

Inis  completes  the  oescript.on  of  the  MAIN  program  except  for  various  print 
flags. 


aaaaMsaaia 


3.0 


SUBROUTINE  ^ROGW  DESCRIPTION  AND  MATHEMATICS 


r 

3.1  Subroutine  MV 

This  is  a double  precision  subroutine  which  premul ti pi ies  a 3 x 1 vector 
by  a 3 X 3 matrix,  "he  ca'l  ’s  MV(M,V,0)  where  M is  a 3 x 3 matrix  with 
real -valued,  double  pr.icisioii  entries,  V and  0 are  3 x 1 vectors  with  real- 
valued,  double-prec's  ion  entries.  M and  V are  the  input,  while  0 is  the 
output. 

3.2  Subroutine  TMEE 

This  is  a double  precision  subroutine  which  converts  elapsed  time  (seconds 
from  0.0  hour,  days,  hour,  minutes,  and  decimal  minutes  in  both  UTC 

and  local  time.  This  subroutine  is  called  at  output  time  only  as  all 
computations  are  referenced  to  UTC.  The  equations  are: 

LD  = UTC  day 

GH  = UTC  hour  = l.T/3600. , where  LT  is  the  time  since  0.0  hour,  UTC, 
in  seconds 

GH  is  checked  to  see  that  it  is  not  less  than  0.  nor  greater  than  or 
equal  to  2T.  If  so,  LD  is  increased  or  decreased  by  1 and  GH  is 
decreased  or  increased  by  24. 

IHG  = integral  part  of  GH 

MG  = (GH  - IHG)  * 60.  = minutes  and  decimal  minutes,  UTC 

Now  local  hour  is  computed.  InitiaT'y,  IDL  = LD,  where  IDL  is  the  loca 

day.  Previously,  the  difference  between  "local"  time  and  UTC  time  for  each 

t h 

site  was  computed,  for  the  K sHe  (the  site  for  which  look  angles  are 
being  computed), 

ILOC  = integral  part  of  (DL0NG(K)-7 . 5)/l 5 
where  DLONG(K)  is  the  W^st  longitude  of  the  site. 

ILOC  effectively  divirP's  tht*  Earth  into  24  time  /ones  with  meridian  of 
Greenwich  cent<>rccl  m /otic*  0. 
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Then, 


IF  ZON  ^ 0. , ILOC  = ILOC  - 24 

This  equation  accounts  for  the  International  Date  Line  and  all  zones  earlier 
than  UTC  (i.e.,  Greenwich). 

Finally, 


TLOC(K)  = ILOC  + TLC 


where 


TLC  = "Legal"  local  time  - "True"  local  time 

See  Figure  1 for  the  "legal"  local  time  zones  on  the  Earth,  excluding  "daylight 
savings."  Note  that  the  TLC  (as  read  in)  can  accomodate  any  time  zone  as  it  can 
be  positive  or  negative  and  has  at  least  three  decimal  digits.  (See  operating 
instructions  in  Appendix  I.) 

Now,  the  local  hour  (XLH)  is  obtained  by: 

XLH  = GH  - TLOC(K) 

XLh  is  then  checked  for  being  less  than  0.  or  greater  than  or  equal  to  24. 

If  true,  IDL  and  XLH  are  changed  in  a manner  similar  to  the  changes  in  LD 
and  GH. 


b 


Following  the  checks  discussed  above. 
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IHL  = integral  part  of  XLH 
ML  = (XLH  - IHL)  * 60 

= Minute  and  decimal  minutes  of  local  time  (ML  is  a double 
precision,  real-valued  variable.) 

The  return  statement  is  then  reached.  The  call  statement  is: 

TIMEE(LD,LS,K,AMG,AML). 

where  the  input  is: 

LD  = day  of  look  angle  epoch, 

LS  = elapsed  time  in  seconds  from  0.0  hour  of  LD,  and 
K = site  index  (K^*’  site). 

The  output  is: 

AMG  = minute  and  decimal  minutes  of  UTC  time,  and 
AML  = minutes  and  decimal  minutes  of  local  time. 

1 

Also,  the  labeled  common  TIME  carries  IDG,  IDL,  IHG,  IHL,  MG,  ML,  and  TLOC(K) 

I from  main  to  the  subroutine  and  back. 

j 

i 3.3  Subroutine  Vernal 

I Ti.is  subroutine  calculates  the  East  longitude  of  the  mean  Vernal  Equinox  at 

i 0.0  hour,  UTl , on  the  given  cay.  This  subroutine  is  taken  from  Satellite 

I liumi nation  Prediction  Using  NORAD  2-Card  Element  Sets,  PAl 00-74-22.  The 

i equations  are: 

1 

! 

I LPYR  integral  part  of  (lYR  - 1973)/4,  where  lYR  is  the  current  year 

U = 366.  * (lYR  - 1973.)  + LPYR  + DAY,  where  DAY  is  the  current  day 


number 


D1  = D + 3400. 


ELNG  = ?60.48673053D0  - 9.85647350  * lO'^  * D - 2.9015  * 10"^^  * 01  * 01 

ELNG  is  then  converted  to  an  angle  less  than  360°  and  greater  than  or  equal  to 
0°.  Then  ELNG  is  converted  to  radian  measure  and  returned. 

ELNG  is  the  East  longitude  of  the  Vernal  Equinox  (sometimes  called  Point  of 
Aries.  In  any  case,  it  is  the  direction  of  the  positive  inertial  X-axis.) 
at  O.U  hour,  UTl , whose  day  number  is  equal  to  DAY. 

The  call  is  VERNAL( lYR, DAY .ELNG) . The  inputs  are  DAY  and  lYR,  while  the 
output  is  ELNG. 

3.4  Subroutine  Block  Data 

This  subroutine  initializes  the  named  commons  of  TIME,  CONST,  YR,  and  INPUT. 

It  is  a double  precision  subroutine. 

The  common  INPUT  consists  of  the  variable  AI,  which  has  dimension  369  x 5. 

It  contains  the  Bessel ian  day  numbers  A,  B,  C,  D and  the  Independent  Day 
number  f as  contained  in  The  American  Ephemeris  and  Nautical  Almanac  for  this 
year.  These  constants  are  referenced  to  the  equinox  of  the  beginning  of  the 
year  until  July  1,  when  the  reference  becomes  the  equinox  at  the  beginning  of 
next  year.  (This  "year"  is  actually  the  "Bessel ian  year",  rather  than  the 
"calendar  year.")  So,  in  addition  to  day  numbers  for  every  day  of  the  year, 
there  are  day  numbers  for  January  0,  December  32,  and  two  sets  for  July  1. 

Except  during  leap  year,  the  369^*^  entries  are  zero  for  all  five  variables. 

All  of  these  numbers  are  supplied  by  data  statements  in  this  subroutine. 

The  common  TIME  has  the  variables  IDG,  IDL,  IH6,  IHL,  XM6,  XML,  TIM,  and  TLOC. 

TLOC  has  dimension  10.  TIM  is  described  in  the  Section  on  Main,  while  the 
others  are  described  in  the  secion  on  the  subroutine  TIMEE.  Block  data  is 
used  to  initialize  IDG,  IDL,  IHG,  IHL,  XMG,  and  XML  to  0. 

The  coimxjn  YR  has  the  data  for  the  current  year.  The  variables  are  XMJYR,  lYR, 
NYR,  and  July.  XMJYR  is  the  modified  Julian  date  of  January  0 of  the  current 
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year  (December  31  of  Iasi.  year).  lYK  is  tin*  year  numbfr  of  the  (urrent.  year, 
and  NYR  is  the  number  of  days  in  Ifiis  yc.ir  . ilUI  Y is  LIk-  day  number  (day)  of 
July  1 for  the  current  year.  All  of  these  variables  are  initialized  via  data 
statements  in  BLOCK  DATA. 

The  last  common,  CONST,  contains  must  of  the  Larth  model  (WGS-72),  time,  and 
numerical  constants  used  in  this  program.  They  are; 

PI  - the  constant  -n  (^3.  IT  159266359873200) , 

RAD  - the  number  of  degrees  in  a radian  (-57. 29377951 3082D0) , 

FILM  - the  number  of  ineters  in  u,i  International  foot  ( = 3.0480-1  ), 

EC2  - the  eccentricity  squared  of  the  Earth  through  the  poles 
(=6.6943177780-3), 

ES  - the  mean  current  eccentricity  of  the  Earth's  orbit 
(=1.67260-2), 

FLAT  - the  flattening  constant  of  the  Earth  (=298.2600), 

AE  - the  average  radius  oi  the  Earth  in  meters  (=637813500), 

WE  - the  average  spin  velocity  of  the  Eartii  with  respect  to 
the  (moving)  mean  Vernal  Equinox  in  radians/second 
(=7.292115854680-5)  for  epoch  1900.0, 

CABER  - the  aberrational  constant  of  the  Sun  (=20.49600),  and 

EPH  - the  current  differer.„e  between  2ph.*meris  time  and  UTC 
(ET  - UTC)  in  se'.onc  > (tPH  car.  oe  found  on  page  >77 
of  the  American  Ephomeris  and  Noutical  Almanac  for  the 
,,ear  desired.) 

3.5  Subroutine  A5>2C.1K 


This  is  a double  precision  subroutine  Wrich  converts  Right  Ascension  and 

1 

I declination  to  X,  Y,  Z direction  cosines  and  vice  versa.  The  flag  input 

is  N ;f  ?'  = 1,  X-Y-7  cosines  re  converted  to  Right  Ascension  and  declination. 
I I-^  N j,  then  Right  Asiension  and  uec  > ii.atiori  are  converted  to  X-Y-Z  cosines. 


i 


Otherwise,  RAS  and  DLC  are  the  inputs  while  X is  the  output.  N is  always 
an  input. 


3.6  Subroutine  CNSTNT 

This  is  a double  precision  subroutine  which  computes  the  various  constants 
needed  in  the  procjrain.  Those  constants  are: 

E - Earth's  mean  obliquity  (dihedral  angle  between  the  equatorial 
plane  and  Earth's  orbital  plane  - called  the  ecliptic  plane). 
This  parameter  is  always  referenced  to  the  beginning  of  the 
nearest  Bessel  ■’an  year  via  use  of  the  flag  LQ. 

XMU  unnuol  genera^  precession  in  declination,  referenced  to  the 
beginning  of  the  Bessel i an  year  (similar  to  E). 

The  equations  for  these  parameters  are: 

TAU  = lYR  - 1900. 

IF(LQ.EQ.O)  GOTO  3 

UB  --  15019.81352  + (365.24219879  - 8.56  * lD-9  TAU)  * TAU  - XMJYR 

DL  = 365.2422  - 1.480-3  * TAU/8.64D-4 

TTO  - iXMJvR  . 1501 9. 499995)73652500 

DAY  = .'’''IJYR  - 15019.813‘'2  -*  OB 

IF'J.GE.JIJLV)  day  - OL  ^ DAY 

DAYl  = DAY  + 3.13520-1  ' 5.0-6 

TT  = DAY/36525. 

TET  = DAYl/36525. 

E = (23.A522944  - 1.301250-2  * TET  - 1.6390-6  * TET  * TET 
I 5.0280-7  * (TFT  ■‘'*3))/RAD 


20 


3 CONTINUE 


1 The  flag  LQ  is  reset  to  LQ  = 1 for  each  run.  Above, 

'I 

j 

j DB  = days  and  decimal  days  from  0.0  hour,  ET,  January  0,  of  the 

current  year  to  the  beginning  of  the  current  Besselian  year, 

OL  = length  of  current  Besselian  year, 

■j 

1 XMJYR  = Modified  Julian  date  of  January  0 for  the  current  year, 

; TTO  = fraction  of  a ccmputat ional  tropical  century  which  has  elapsed 

since  the  beginning  of  the  Ephermeris  Julian  year  of  1900, 

TT  = fraction  of  a computational  tropical  century  which  has 

elapsed  since  the  beginning  of  the  Ephemeris  Besselian  year 
of  1900, 

TET  = TO, 

UTC,  E,  XMO,  and  LQ  were  previously  defined. 

Now  the  Besselian  constants  for  me  current  day  (day  number  J)  are  taken  from 

Block  Data  and  irans.  erred  to  the  variables  AA,  uA,  and  O'^A  --  all  of  dimension 

5.  OA  and  D2A  contain  the  first  and  se< ond  differences,  respectively,  of  the 

uesselian  constants  A o’'  these  variables  are  in  radians,  radians/day,  or 
2 

r .diai.s/day  , as  i 'Cc  it  equatiens  are: 

J1  = J + , 

IF  (J.GE.  Jui'/l  J!  “ Ji  + I 
AA(I)  - A..^.;l  ..r)/ractor;  1 = 1,0 

Now  1 "*  rep.'S.  't  t?  c.  ' nt  lay  yesterday,  and  "1"  tomorrow, 

en: 

DA(0)  = y-''  .)  - AA(0; 

= .el)  + o\(-l)  -2  * AA(0)  for  the  five  AA(  )'s 
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Also,  AA(I),  and  D2A(I)  are  calculated  for  day  J + 1 and  named  AAP, 
and  D2AP,  respectively,  but  DAP  = DA. 


• • 


These  constants  are  used  in  the  following  manner.  Let  XL  = fraction  of  day  J 
for  which  look  angles  are  calculated.  Then: 

C(I)  = AA(I)  + DA(I)  * XL  + D2A(I)  * XL  * XL/2.,  I = 1,5 

if  XL  is  less  than  or  equal  to  1/2.  Otherwise, 

C(I)  = AAP(I)  - DAP(I)  * (1  - XL)  + 

D2AP(I)  * (1  - XL)  * (1  - XL)/2.,  I = 1,5 

This  computation  is  accomplished  in  the  MAIN  program. 

Due  to  the  discontinuity  at  July  1 and  the  end  of  the  year,  the  AAP,  DAP,  and 
D2AP  vectors  are  reset  for  those  two  days. 

IF  ((J.NE.(July  - 1).AND.{J.NE.NYR))  GOTO  2 
AAP(I)  = AA(I) 

DAP(I)  = -DA(I) 

D2AP(I)  = D2A(I) 

2 CONTINUE. 

Also,  in  MAIN,  XL  is  used  past  XL  = 1/2  instead  of  using  1 - XL. 

Finally,  tne  parameters  EOFE,  DEFE,  and  DW  are  calculated.  These  are: 

EOFE  = "Equation  of  Equinoxes"  (difference  between  mean  and  apparent 
Vernal  Equinox  at  0.0  hour,  UTC), 

OFFE  = difference  in  EOFE  from  day  J to  day  J+1 , and 

DW  = secular  change  in  Earth  spin  rate. 
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The  equations  are: 


TT1  = TTO  ^ J/36525D0 

TTAU  = (J  - DB)/365.2422UU 

TTAUI  = TTAU  + 1 . /365. 2422D0 

IF  (J.GE.JULY)  TTAU  = TTAU  - t)L/365.2422DO 

IF  (J.GE.JULY)  TTAUI  = TTAUl  - DL/365. 242200 

XM  = XMO  * TTAU 

XMl  = XMO  * TTAUI 

EOFE  = (AA(1)  - XM)/ETN 

EOFEl  = (AAP(l)  - XM1)/ETN 

DEFE  = EOFEl  - EOFE 

DW  = 4.28  * 10-15  * TTl 

Note  that  AA(1)  is  the  Bessel i an  Day  Number  A. 

The  call  to  CNSTNT  is  CNSTNT(J ,E ,OW,EOFE ,OEFE ) . 

The  only  input  here  is  J,  the  day  number.  Other  inputs  obtained  from 
common  are: 

Coiimon  CONST  - the  constant  RAO, 

Conmon  YR  - the  yearly  constants. 

Common  INPUT  - the  Bessel ian  Day  Numbers,  and 
Coiimon  (unnaiiied)  - the  variables  FACTOR,  OSNl  , and  LQ. 

FACTOR  is  the  number  of  arcseconds  iri  a radian,  while  05N1  is  the  double 
precision  sine  of  an  arcsecond. 

Outputs  are  E,  DW,  EOFE,  and  DEFE.  Also,  the  current  Besselian  Day  Numbers 
and  those  for  the  succeeding  day,  together  with  their  first  and  second 
differences,  are  passed  through  the  named  cotmon  NOW.  The  fraction  of  a 
tropical  year  which  rtas  elapsed  (will  elapse)  since  the  beginning  of  the 
nearest  Besselian  year,  TTAU,  is  passed  through  the  unnamed  common. 
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There  is  a second  entry  to  CNSTNT,  called  CNSTJ.  The  call  is 
CNSTJ(J,IJ,RAS,DEC,XJJ,XJP).  The  inputs  are: 

J - the  current  day  number; 

IJ  - set  in  MAIN,  IJ  = -1  if  the  star  declination  is  negative, 

IJ  = +1  otherwise; 

RAS  - updated  mean  Right  Ascension  of  star;  and 
DEC  - updated  mean  declination  of  star. 


The  outputs  (in  units  of  radians)  are: 

XJJ  - the  2^^  order  correction  to  Right  Ascension,  and 
XJP  - the  2^*^  order  correction  to  declination. 


The  equations  are: 

PI  = (A(J)+IJ*D(J))*DSIN{RAS)+(B(J)+IJ*C(J))*DCOS(RAS) 
P2  = (A(J)+IJ*D(J))*[)COS(RAS)-(B(J)+IJ*C(J))*DSIN{RAS) 
XJJ  = PI *P2*DSN1 /FACTOR 
XJP  = -Pl*Pl*DSNl/(2. ♦FACTOR) 


Above,  A(J),  B(J),  C(J),  and  0(J)  are  the  Besselian  day  numbers  A,  B,  C,  and 
0 for  day  J. 


3.7  Subroutine  OCT 


This  subroutine  changes  degrees  into  octal  encoder  readout.  Some  explanation 
of  encoder  readout  is  required.  The  encoder  has  a set  number  of  "bits"  - 
say,  17.  Each  bit  is  binary  (either  on  or  off).  The  first  bit  (reading 
from  left  to  right)  is  "1"  if  the  angle  is  above  180®.  Otherwise,  it  is 
"0".  The  second  bit  is  "1"  if  the  angle  (or  the  remainder  of  angle  -180”  if 
bit  1 = 1)  is  above  90®  (90  = 180/2),  otherwise  it  is  "0".  The  third  bit 
is  "1"  if  the  angle  (for  the  appropriate  remainder)  is  above  45®  (45  * 90/2) 
and  is  "0"  otherwise.  This  process  continues  down  to  the  least  significant 
bit  (LSB).  In  this  case,  the  17th  bit  = 9.89  arcseconds.  See  Table  lA  for 
a short  summary  of  LSB  values  for  popular  bit  configurations. 
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TABLE  lA 

Short  Table  of  LSB's  for  Common  Angle  Encoder  Bit  Configurations 


BIT  # 

LSB 

ARCSECONDS 

MILS 

13 

158.202 

0.7825, 

15 

39.548 

0,19531 

17 

9.889 

0.04883 

19 

2.473 

0.01221 

21 

0.619 

0.00305 

23 

0.155 

0.00076 

After  determining  the  binary  bit  number,  every  set 

of  three  starting  with 

the  LSB  is  changed  to  octal.  As  an  example,  let  us  change  the  following 
angle  to  17  bit  octal  representation. 

TABLE  IB 

Computation  of  Binary  17-B1t  Representation 

Angle  = 123.426“ 

BIT  # 

VALUE 

REMAINDER 

1 

0 

123.426 

2 

• 1 

43.426 

3 

0 

43.426 

4 

1 

20.926 

5 

1 

9.676 

6 

1 

4.051 

7 

1 

1 . 2385 

8 

0 

1 . 2385 

9 

1 

0.535375 

10 

1 

0.183812 

11 

1 

0.008031 

12 

0 

0.008031 

13 

0 

0.008031 

14 

0 

0.008031 

15 

0 

0.008031 

i6 

1 

0.002538 

17 

0 

0.002538 

A 


— . ' ' ' "1^  BP 

The  binary  value,  in  groups  of  3,  is:  01  011  110  111  000  010.  The  octal 
representation  of  that  number  is  136702. 

The  process  as  shown  above  is  similar  to  the  process  used  in  OCT.  First,  there 
is  a second  entry  called  OCTO.  The  call  is  OCTO(IBITE).  The  input  is  IBITE, 
the  number  of  encoder  bits,  read  from  the  station  card.  If  IBITE  = 0,  this  call 
is  skipped.  In  OCTO  the  number  of  groups  of  three  are  determined,  as  are  the 
number  of  bits  in  the  first  set  (3,  2,  or  1).  The  division  for  the  first  group 
of  three  is  thus  determined  (180°  if  the  number  in  the  first  set  is  3,  90°  if 
it  is  1,  and  45°  if  it  is  2).  Several  other  flags  are  also  set  to  be  used 
when  OCT  is  cal  led. 

The  call  for  OCT  is  0CT(X,I0UT),  where  X is  the  input  - a double  precision 
number  in  decimal  degrees  - and  lOUT  is  the  octal  (integer)  representation. 

The  heart  of  OCT  is  the  double  DO  loop  shown  below. 

DO  6 K = LO,  I 
DO  5 J = 1 ,3 
R = 2.**J 

IF(XF.GE.(DIV/R))  IBIT(K)  = IBIT(K)  + (8./R) 

5 IF(XF.GE.(D1V/R))  XF  = XF  - (DIV/R) 

6 DIV  = DIV/8D0 

Above,  XF  is  the  remainder  passed  into  this  loop  after  processing  the  odd  number 
of  front-end  bits.  DIV,  LO,  and  L are  set  in  OCTO.  The  vector  IBIT  stores  the 
result.  IBIT  is  zeroed  on  entry  into  OCT. 

After  this  processing  is  complete,  the  remainder  is  checked  to  see  if  it  is 
greater  or  equal  to  three-fourths  the  LSB.  If  so,  the  last  octal  number  is 
increased  by  1.  (The  number  "3/4"  is  an  "engineering"  choice  based  on  the 
inherent  noise  in  the  LSB  of  most  encoders.  A mathematical  choice  would  be 
"1/2".) 
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If  the  smallest  octal  Integer  Is  Increased,  It  Is  checked  to  see  if  it  is  8. 
If  so.  it  is  zeroed  and  the  next  smallest  integer  is  increased  by  one.  This 
continues  until  an  integer  is  reached  which  is  not  8 after  increase  or  until 
the  first  integer  is  increased.  It  is  then  checked  to  see  if  the  number  is 
equivalent  to  zero  in  octal  representation.  If  so,  the  entire  number  is 
zeroed.  Figure  2 shows  an  example  of  many  types  of  octal  output.  Figure  3 
shows  four  examples  in  the  Azimuth  Octal  column  of  the  roundup  portion  of  the 
subroutine. 

Finally,  the  single-digit  integers  in  the  vector  IBIT  are  converted  into  a 
multi-digit  integer.  This  number  is  lOUT.  IF  X,  the  input,  is  negative, 
then  lOUT  is  multiplied  by  -1. 
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4.0 


EFFECT  OF  CORRECTIONS  AND  APPROXIMATIONS 


In  this  section  the  effect  of  the  various  corrections  to  the  mean  star 
position  (transformation  from  "mean"  to  "apparent")  and  to  the  apparent 
star  position  (transformation  from  "apparent"  to  "look  angles")  will  be 
discussed.  This  section  consists  of  self  comparisons  of  the  verified  program. 
Its  purpose  is  to  establish  a basis  for  estimating  the  numerical  "size" 
of  various  corrections  before  discussing  verification  efforts.  These  com- 
parison were  made  by  use  of  the  "degrees"  output,  which  has  a LSB  of  0.00001 
degrees  or  0.036  arcseconds.  Hence,  the  tables  in  this  report  only  show 
changes  of  about  0.04  arcseconds  or  more.  This  step  size  is  only  in  the 
printout  and  not  in  the  computer  processing. 


4.1 


Transfonnation  from  "MEAN"  to  "APPARENT" 


The  transformations  here,  in  order  of  use,  are: 


a.  Mean  motion  is  corrected  for  "proper  motion''  with  respect 


to  the  mean  equator  and  equinox  of  the  nearest  Bessel ian  New 
Year; 

b.  The  rectangular  coordinates  are  corrected  for  "annual" 
aberration; 

c.  The  rectangular  coordinates  are  transformed  to  the  "true" 
equinox  and  "true"  equator  of  date  (precession-nutation  corrections); 

d.  The  annual  parallax  correction  (transformation  from  heliocen  - 
trie  to  geodetic  coordinates)  is  introduced;  and 

e.  Second  order  corrections  caused  by  cross-terms  involving 
aberration  corrections  for  a fixed  star  system  and  precession  - 
nutation  corrections  for  a moving  star  system  are  introduced. 

The  assumptions  involved  are: 

a.  No  star  has  a declination  of  exactly  + 90°. 

b.  Second  order  proper  motion  corrections  are  ignored. 

This  is  tenable  because  the  largest  corrections  used  here 
cover  only  1/2  year,  while  second  order  proper  motion  correc- 
tions are  on  the  order  of  0.01  arcsecond/tropical  century. 

in 


1 


I 


,! 

} 


c.  "&n«ll  angle"  approximation  (s1n(X)  - X and  cos(X)  ■ l.-X^/2., 
where  X Is  In  radians  and  all  second  order  terms  are  omitted  In 

the  final  computation)  is  used  throughout  this  program,  particularly 
in  the  precession-nutation  computation.  Replacement  of  these 
approximations  by  actual  sines  and  cosines  resulted  in  no  difference 
to  the  order  desired  (5x10  arcseconds). 

d.  The  parallax  correction  is  omitted  if  tt  < 0.01  arcseconds. 

Also,  this  correction  is  assumed  constant  over  a day's  period. 

Again,  the  error  is  less  than  5x10  arcseconds. 

e.  One  can  also  introduce  second  order  errors  for  the  precession- 

f 

nutation  transformation  (see  reference  1,  p.  370).  However,  there 
is  no  evidence  that  this  is  required  for  0.01  arcsecond  accuracy. 

f.  Corrections  for  the  Earth's  elliptic  terms  are  not  used.  These  . 

corrections  amount  to  0.001  arcsecond/tropical  century  and  are  not  ? 

significant.  \ 

g.  Additional  assumptions  are:  1 

I 

1.  The  motion  of  the  Earth  lies  in  the  ecliptic;  I 

2.  Second  order  effects  of  aberration  are  negligible;  and  1 

3.  "Light"  time  and  "secular"  aberration  are  never  corrected 
(for  stars). 

The  following  effects  are  arranged  in  order  of  decreasing  effect.  The  process 
followed  is  to  compare  outputs  run  from  a deck  with  all  effects  "zeroed  out" 
against  a program  run  with  all  effects  "zeroed  out"  except  the  stated  effect. 

The  effects  are: 


Diurnal  aberration; 

Update  of  Bessel ian  constants  by  first  and  second  differences; 

of  Bessel  ian  constants  about  a period  of  [0,1]  instead 
of  o period  of  [-1/2,  1/2]  - the  fractions  refer  to  "days"; 
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Update  of  Bessellan  constants  using  first  differences  only; 

Linear  update  of  Equation  of  Equinoxes; 

Second  order  calculations  (J  and  J'  numbers);  and 
Secular  change  to  Earth  spin  rate. 

4.2  Diurnal  Aberration 

The  order  of  all  the  tables  will  be  "corrected  - uncorrected." 

Table  3 below  Is  "Diurn  Corrected"  - "Std."  Table  2 lists  the  Input 
data,  which  Is  the  same  for  the  following  tables,  except  for  day  of  run. 

For  Table  3 the  "rough"  azimuth  and  elevation  for  the  "std"  run  are  listed 
for  help  In  Interpreting  the  data.  The  other  tables  for  this  same  day  would 
have  the  same  "rough"  pointing  angles. 

TABLE  2 

Input  Data  for  Runs  Used  In  Section  4 

Latitude  (site)  = 37.4958°  N 
Longitude  (site)  = 237.4961°  E 
Height  above  ellipsoid  = 0.  ft. 

Astronomic  latitude  - Geodetic  latitude  = 0. 

Astronomic  longitude  - Geodetic  longitude  = 0. 

Minimum  printed  elevation  = 0.° 

Star  Name  = ALFA  BETA^ 

Star  Mean  declination  = 56°24'22" 

Star  Mean  Right  Ascension  = 99°47'1.5"(6  hours, 39  m1n,8.1  sec) 

Star  Epoch  = 1976.0 

Month  and  Day  of  Run  =1/1  GMT 

Note  that  the  DIURN  parameter  (In  arcseconds)  Is: 

DIURN  = 3.198D-1*(AEL+HT)*DC0S(DLAT), 

where  (AEL  + HT)  Is  the  distance  (In  standard  Earth  radii)  from  the  site 
to  the  Earth's  center  and  DLAT  Is  the  geodetic  latitude.  This  parameter 
has  a maximum  of  0.32  arcseconds  for  HT  ■ 0 ft.  and  DLAT  • 0°. 

t Not  listed  In  most  star  catalogs. 
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TABLE  3 

Diurnal  Aberration  Test 


4.3 


Update  of  Bessel i an  Constants 


The  second  effect  studied  is  the  use  of  first  and  second  differences  to 
update  the  Bessel ian  constants.  In  Table  4 the  differences  between  "usual"  - 
"no  update"  are  recorded.  Note  that  these  differences  are  large  enough  to 
require  interpolation  for  high  accuracy  work. 


j 


TABLE  4 

Test  for  Full  Updating  of  Bessel ian  Constants 


Time 

GMT 


Azimuth  Diff.  Elevation  Diff. 

arcseconds  arcseconds 


1:00-1/1 

0. 

2:00 

-.04 

3:00 

-.04 

4:00 

-.04 

5:00 

-.07 

6:00 

-.04 

7:00 

+.04 

8:00 

+ .18 

9:00 

+ .29 

10:00 

+ .29 

11:00 

+ .22 

12:00 

+ .14 

13:00 

+ .11 

14:00 

+ .04 

15:00 

+ .04 

16:00 

0. 

17:00 

-.04 

18:00 

-.07 

19:00 

-.14 

20:00 

-.18 

21:00 

-.22 

22:00 

-.25 

23:00 

-.29 

0:00-1/2 

0. 

1:00 

0. 

0. 

0. 

-.04 
-.04 
-.04 
-.07 
-.07 
-.11 
-.04 
+ .04 
+ .04 
+ .11 
+ .14 
+ .18 
+ .18 
+ .18 
+ .1C 
+.22 
+ .18 
+.18 
+.14 
+ .11 
+ .04 
0. 

0. 
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Below,  Table  5 records  the  differences  between  the  look-angles  produced  by 
the  usual  update  method  (update  Interval  being  [-1/2,  + 1/2)  in  terms  of 
days  about  the  midnight  of  the  interpolated  day)  to  those  produced  by 
use  of  an  interpolation  period  of  [0,1).  The  differences  are  "[0,1)  period"  - 
"usual".  Table  6 records  the  differences  between  the  usual  method  (first  and 
second  differences)  and  the  values  produced  by  use  of  first  differences  only. 
Again,  the  differences  are  recorded  in  the  sense  "first  only"  - "usual." 


TABLE  b 

Test  for  Updating  of  Bessel ian  Constants, 
Change  of  Usual  Interpolation  Interval 


Time 

GMT 

1:00-1/1 

2:00 

3:00 

4:00 

5:00 

6:00 

7:00 

8:00 

9:00 

10:00 

11:00 

12:00 

13:00 

14:00 

15:00 

16:00 

17:00 

18:00 

19;0C 

■; 

L I : )'j 
22:00 
2a:00 
j:C  ■./2 

1 :or. 


Azimuth  Diff. 
arcseconds 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

+ .04 
+ .04 
+ .11 
+ . 14 
+ .18 
+ .22 
+ .25 
0. 

0. 


Elevation  Diff. 
arcseconds 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-.04 

-.07 

-.07 

-.07 

-.11 

-.14 

-.14 

-.14 

-.11 

-.07 

-.04 

0. 

0. 


3b 


TABLE  6 

Test  for  Updating  of  Bessel ian  Constants, 
1st  Differences  Only 


Time 

GMT 


Azimuth  Diff.  Elevation  Diff. 

arcseconds  arcseconds 


1:00-1/1 

0. 

2:00 

0. 

3:00 

0. 

4:00 

0. 

5:00 

0. 

6:00 

0. 

7:00 

0. 

8:00 

0. 

9:00 

0. 

10:00 

0. 

11:00 

0. 

12:00 

0. 

13:00 

0. 

14:00 

0. 

15:00 

0. 

16:00 

0. 

17:00 

0. 

18:00 

0. 

19:00 

0. 

20:00 

0. 

21:00 

0. 

22:00 

0. 

23:00 

0. 

0:00-1/2 

0. 

1 :00 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

+ .04 

0. 

0. 

-.04 

-.04 

0. 

0. 

0. 


0. 

0. 


0. 

0. 

0. 


0. 

0. 

0. 

0. 


From  Table  5 one  can  observe  that  the  differences  increase  with  time. 

Also,  in  the  original  runs  the  differences  decreased  with  smaller 
elevation  angle.  Note  that  these  differences  are  much  larger  than  those 
caused  by  omitting  the  second  order  interpolation.  The  differences  caused 
by  use  of  second  order  interpolation  as  seen  in  Table  6 are  small.  These 
differences  will  rarely  be  significant.  See  the  graph  of  second  order 
coefficient  values  in  the  next  section. 


4.4 


Small  Effects 


A.  Linear  Change  to  Equation  of  Equinoxes 

The  parameter  known  as  Equation  of  Equinoxes  is  the  difference  between  the 
mean  longitude  of  the  Vernal  Equinox  at  0:00  hour,  UTl , and  the  "True" 
longitude  (i.e.,  inclusion  of  periodic  terms).  This  difference  is  on  the 
order  of  0.8  seconds  and  changes  by  about  .005  seconds  per  day.  Since  the 
updated  Bessel ian  constants  transform  the  coordinates  to  the  instantaneous 
"true"  equinox,  the  Earth  spin  must  be  referenced  to  that  same  equinox. 
This  is  achieved  by  the  following  equations: 

E0FE2  = EOFE  + DEFE  * XL11 
ELNG  = ELNGO  + E0FE2 

where 

ELNG  = "True"  longitude  of  the  Vernal  Equinox, 

ELNGO  = "Mean"  longitude  of  the  Vernal  Equinox, 

EOFE  = "True"  - "Mean"  longitude  at  the  same  epoch, 

DEFE  = first  difference  of  EOFE  (per  day), 

XLll  = fraction  of  elapsed  day  since  0:00  hour  (UTl  time),  and 
E0FE2  = updated  Equation  of  Equinoxes. 

The  differences  between  Equation  of  Equinox  "Corrected"  minus  "Standard" 
(no  corrections)  as  reflected  in  the  look  angles  are  listed  in  Table  7. 

,\  hand  check  of  the  2nd  differences  for  the  Equation  of  Equinoxes  showed 
k„-,ac  their  siiie  aid  not  warrant  their  use  in  this  program. 
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TABLE  7 

Test  for  Linear  Update  to  Equation  of  Equinoxes 


Time  Azimuth  Diff. 

GMT  arcseconds 


1:00-1/1  0. 

2:00  0. 

3:00  +.04 

4:00  +.04 

5:00  -.04 

6:00  -.04 

7:00  -.04 

8:00  -.07 

9:00  -.04 

10:00  -.04 

11:00  0. 

12:00  0. 

13:00  0. 

14:00  0. 

15:00  +.04 

16:00  +.04 

17:00  +.04 

18:00  +.04 

19:00  +.04 

20:00  +.07 

21:00  +.07 

22:00  +.04 

23:00  +.04 

0:00-1/2  0. 

1:00  0. 


Elevation  Diff. 
arcseconds 

0. 

0. 

0. 

0. 

+.04 

0. 

0. 

0. 

0. 

-.04 

-.04 

0. 

-.04 

-.04 

-.04 

-.04 

-.04 

-.04 

0. 

0. 

0. 

+ .04 
+ .04 
0. 

0. 
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B.  Earth  Spin  Rate  Changes 


The  second  small  change  is  In  the  Earth's  spin  rates.  The  constant  part, 
WE,  Is  the  value  for  epoch  1900.  There  is  a secular  increase,  DW,  whose 
value  is: 

T = time  in  days  since  1900 
TT  = T/36525 

DU  = 4.29  D-15  * TT  in  radians/second 

Then  the  corrected  spin  rate  (relative  to  the  moving  Vernal  Equinox) 
is: 

WEOl  = WE  + DU 

Here,  WE  = 7.29211585468  * 10"^  radians /second,  and 
DU  ^ 3.25  * lO"^'*  for  1976. 

However,  there  is  another  correction  due  to  polar  motion.  This  is  a scale 
factor  of  the  form: 

WED  = WEDl  * (1.  + DUTDOT) 

where  DUTDOT  is  the  rate  of  change  of  UTl  in  the  dimensionless  form  of 
days/day.  DUTDOT  is  preset  to  -2.0  mil  1 iseconds/day  (unless  changed  by 
polar  motion  corrections).  Hence, 

DUTDOT  = -2.0  mil  1 iseconds/day  = -2.3  * 10'®  days/day,  and 
i.  + DUTDOT  = .999999977. 

I";  p-.  ..  ,cu .c:r,  WED  - WEDl  = -2.79  * 10'^^  radians/second.  Hence,  the 
p .ridpal  erfort  is  that  of  the  motion  of  the  spin  axis  (polar  motion 
" %?)  at  the  present  time.  The  secular  change  is,  or  course,  progressive. 
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effect  of  this  on  look  angles  amounts  to  no  more  than  0.04  arrseconds 


dur  iric 

any  day.  See  Table 

' 8 for  a list 

of  these  differences 

("corrected"  - 

"uricoi- 

reefed",  as  usual). 

TABLE 

8 

Test  for 

Correction  to 

Earth  Spin 

Rate 

r iiiu> 

Az.  Diff. 

El.  Diff. 

Time 

Az.  Diff. 

El.  Diff. 

GMT 

arcsec. 

arcsec. 

GMT 

arcsec. 

arcsec. 

1 .MO- 

1/'  0. 

0. 

13:00 

0. 

0. 

r;00 

0. 

0. 

14:00 

-.04 

+ .04 

3:00 

0. 

0. 

15:00 

0. 

+ .04 

4:00 

0. 

0. 

16:00 

0. 

0. 

r no 

0. 

0. 

17:00 

0. 

0. 

6:00 

0. 

0. 

18:00 

0. 

0. 

7:00 

f.04 

0. 

19:00 

-.04 

0. 

8:00 

+ .04 

0. 

20:00 

0. 

0. 

9:00 

+ .04 

0. 

21:00 

0. 

0. 

10:00 

+ .04 

0. 

22:00 

-.04 

0. 

11:00 

+ .04 

+ .04 

23:00 

-.04 

-.04 

12:00 

0. 

0. 

0:00-1/2 

0. 

0. 

1:00 

0. 

0. 

C.  S * ;V/,i  Order  Effects 

The  second  order  effects  are  the  cross- terms  between  aberration  corrections 
performed  in  a fixed  coordinate  system  and  nutation  effects,  which  occur  in 
a moving  coordinate  system.  These  effects  are  very  small  except  for  stars 

at  high  declination  angles.  In  the  case  considered  above  the  effect  on  azimuth 

is,  at  most,  .04  arcseconds  for,  on  the  average,  one  point  per  day.  However, 
in  Table  9 below,  based  on  a Right  Ascension  of  21  hours  for  day  26,  1976, 

it  is  shown  how  this  effect  grows.  For  stars  with  declination  below  60® 

absolute  value,  second  order  corrections  need  not  be  applied. 
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TABLE  9 


Table  of  2nd  Order  Effects  Versus  Star  Declination 


Declination 

degrees 

Rt.  Asc.  Correction 
arcseconds 

Dec.  Correction 
arcseconds 

0 

0. 

0. 

10 

-.00004 

-.00017 

20 

-.00018 

-.00035 

30 

-.00044 

-.00056 

40 

-.00093 

-.00081 

50 

-.00188 

-.00115 

60 

-.00396 

-.00167 

70 

-.00997 

-.00264 

80 

-.04250 

-.00546 

89 

-4.33669 

-.05513 

-10 

+.00002 

+.00016 

-20 

+.00006 

+.00032 

-30 

+.00016 

+.00051 

o 

1 

+.00035 

+.00074 

-50 

+.00070 

+.00105 

-60 

+.00147 

+.00153 

-70 

+.00370 

+.00243 

-80 

+.01576 

+.00501 

-89 

+1.60842 

+.05061 
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Total  Small  Efforts 


Table  10  below  portrays  the  differences  between  the  uncorrected  data  and  the 
data  with  all  three  small  effects  corrected.  Again,  the  differences  are 
"corrected"  - "uncorrected". 


TABLE  10 

Test  for  Total  Small  Errors  Effects 


I line 

Azimuth  Diff. 

Elevation  Diff, 

Gi4T;_ 

arcseconds 

arcseconds 

1 : 00-1 /I 

0. 

0. 

2:00 

0. 

0. 

3:00 

+ .04 

0. 

4:00 

+ .04 

0. 

b:00 


0. 


+ .04 


6:00 

0. 

7:00 

0. 

8:00 

-.04 

9:00 

-.04 

lo.-oo 

0. 

11:00 

0. 

12:00 

0. 

0. 

0. 

0. 

0. 

0. 

-.04 


0. 


13:00 

0. 

'4:00 

0. 

15:00 

0. 

16:00 

+ .04 

17:00 

+ .04 

18:00 

+ .04 

19:00 

0. 

20:00 

+ .04 

21:00 

+ .04 

22:00 

0. 

23:00 

+ .04 

0. 

0. 

0. 

0. 

-.04 

0. 

0. 

0. 

0. 

+ .04 


0. 


0:00-1/2  0.  0. 

!:'J0  0.  0. 


‘;j 


f 
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4.5  Effect  of  Proper  Motion  and  Parallax 

A.  Parallax 

This  correction  is,  in  general,  small.  In  fact,  it  can  almost  always  be 
neglected  if  unknown.  The  known  values  on  parallax  can  be  found  in  the 
Yale  University  Catalog  and  its  supplements.  A copy  is  in  the  possession 
of  Dr.  George  Sinclair  at  FEC/Performance  Analysis  Department. 

The  largest  parallax  known  has  a value  of  C 762  arcseconds.  However,  only  66 
stars  have  parallaxes  greater  than  or  equal  to  0.100  arcseconds.  Further,  the 
largest  negative  parallax  is  -.057  arcseconds.  In  Table  11,  runs  with  these 
two  parallaxes  are  compared  against  a run  with  no  parallax.  Again,  "difference"  = 
"corrected"  - "uncorrected". 

One  should  note  that,  in  the  expression  for  parallax  correction  (written  below 
in  rectangular  coordinates),  it  is  assumed  that  the  Earth's  orbit  is  circular 
and  unperturbed  and  that  the  Bessel ian  numbers  C and  D are  referenced  to  the 
Sun's  center.  Actually,  the  Bessel ian  numbers  C and  D are  referenced  to  the 
barycenter  of  the  solar  system  (approximately).  However,  these  corrections, 
when  multiplied  by  the  parallax,  are  effectively  of  third  order  and  can  be 
neglected.  The  equations  mentioned  above  are: 

AX  = C*PAR/(CABER*DCOS(E)) 
ay  = D*DCOS(E)*PAR/CABER 
AZ  = AY*DTAN(E) 

where , 

CABER  = 20.496  arcseconds  = solar  aberration  constant, 

- mean  obliquity  of  Earth,  referenced  to  nearest  Bessel ian  New  Year, 


TABLE  ri 

Effect  of  Parallax  Corrections 


P^ral  Ux 

= ^.752 

Parallax 

= -.057 

Parallax  = 

+_J. 

Time  Az.  Liiff. 

arcsecono 

El.  niff, 
^rc^c'^w 

Az.  Diff. 
otT second 

El.  Diff. 
arcsecond 

Az.  Diff. 
arcsecond 

El.  Diff. 
arcsecond 

1 ■ 00- 1 / 1 

-.36 

+ .2C 

-r.04 

0. 

-.04 

+ .04 

2:00 

-.43 

+ 18 

+ .04 

0. 

-.07 

+ .04 

3:00 

-.50 

+ .11 

+ .07 

0. 

-.04 

0. 

4:00 

-.58 

0. 

+.07 

0. 

-.07 

0. 

5:00 

-.72 

-.04 

+ .04 

+ .04 

-.11 

0. 

6:00 

-.79 

-.18 

+ .04 

0. 

-.11 

-.04 

7:00 

-.72 

-.30 

+ .07 

+ .04 

-.11 

-.04 

8:00 

-.14 

-.43 

0. 

0. 

-.04 

-.07 

9:00 

■S-.58 

-.36 

-.04 

+ .04 

+ .07 

-.04 

10:00 

-I-. 79 

-.25 

-.04 

0, 

+ .11 

-.04 

11:00 

-*-.72 

-.14 

-.04 

0. 

+ .•(1 

-.04 

12:00 

+ .61 

0. 

-.07 

0. 

+ .07 

0. 

13:00 

-r.54 

+ . 07 

-.04 

0. 

+ .07 

0. 

14:00 

+ .43 

+ .14 

-.04 

0. 

+ .04 

+ .04 

15:00 

+ .40 

+ .22 

-.04 

0. 

+ .04 

+ .04 

16:00 

+ .30 

+ .29 

-.04 

-.04 

+ .04 

■r,04 

17-00 

+ .22 

+ .30 

0. 

-.04 

+ .04 

+ .04 

18:00 

+ .18 

+.36 

0. 

-.04 

+ .04 

+ .04 

19:00 

+ .11 

+ .40 

-.04 

-.04 

0. 

+ .04 

20:00 

+ .04 

+ .43 

0. 

-.04 

0. 

+ .07 

21:00 

-.07 

+ .43 

0. 

-.04 

0. 

+ .07 

22:00 

-.18 

+ .40 

0. 

0. 

-.04 

+ .07 

23:00 

-.25 

+ .36 

0. 

-.04 

-.04 

+ .04 

0:00-1/2 

-.29 

+ .29 

+ .04 

-.04 

-.04 

+ .04 

1:00 

-.36 

+ .25 

+ .04 

0. 

-.04 

+ .04 

Note  above 

that  the 

corrections 

are  larger  than  the  design 

accuracy  of  0. 

01 

arcseconds . 

Hence, 

if  tlie  parallax  is  known. 

it  should  be 

used  in  STARl 

runs. 
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B. 


Effect  of  Proper  Motion 


'i 

,1 


Proper  motion  changes  the  mean  coordinates  of  the  star  over  a year's  period. 
Hence,  there  will  be  no  effect  from  proper  motion  at  the  beginning  and  end 
of  the  year,  and  there  will  be  a maximum  effect  at  July  1 (mid-year).  Note, 
here,  that  the  reference  epoch  is  changed  at  July  1 so  that  the  size  of  the 
annual  effects  are  minimized.  Table  12  lists  the  differences  ("run  with 
proper  motion"  - "std".)  for  various  combinations  of  proper  motions,  all 
large,  on  day  182.  So,  this  table  can  be  considered  a table  of  "worst"  cases. 
UR  and  UD  are  listed  in  arcseconds/tropical  year. 

Note  the  size  of  the  differences.  In  order  to  preserve  the  accuracy  standards 
of  this  program,  proper  motion  corrections  must  be  included  if  either  UR  or  UD 
is  larger  than  0.1  arcseconds/tropical  year  and  if  the  month  of  the  run  time 
is  between  March  and  October. 


TABLE  12 

Effect  of  Large  Proper  Motion  Corrections 


4.6 


Effect  of  UTl  Correction 


A.  Values  of  Polar  Motion  Coefficients 

The  polar  motion  phenomena  can  be  represented  by  three  numbers.  The  first  Is 
the  difference  in  time  = "UTl"  - "UTC"  (in  seconds).  The  last  two  refer  to  the 
offset  of  the  pole  from  the  CIO  (Conventional  International  Origin)  along  the 
X and  Y axis  of  the  tangent  plane.  These  offsets  are  in  arcseconds.  Sets 
of  these  numbers  are  available  in  the  "US  Naval  Observatory  Bulletin" 

(Series  7),  "BIH  Circular  D119",  and  "Monthly  Notes  of  the  International 
Polar  Motion  Service."  Preliminary  estimates  (before  epoch)  are: 

1.  DUTl  correction.  (Time  only);  accurate  to  +.05  seconds  (with 
respect  to  other  preliminary  estimates).  This  is  changed  monthly. 

2.  Estimate  correction.  This  is  a time  equation,  good  for  60  days 
and  available  from  the  Series  7 bulletin. 

3.  Extrapolated  correction.  Time  only;  this  is  a daily  estimate 
for  a week,  and  is  available  one  week  ahead.  The  rate  of  change 
of  UTl  - UTC  (in  milliseconds/day)  is  also  available  in  this  table 

of  estimates.  | 

1 

Measured  values  available  are:  | 

I 

1.  USNO  Rapid  Service.  Daily  values  of  UTl  - UTC,  X,  and  Y for  | 

a week,  available  one  week  after  the  epoch.  This  is  obtained  from  | 

1 

satellite  tracking.  I 

2.  BIH  (Bureau  International  de  THeure)  values.  Values  of  | 

UTl  - UTC,  X,  and  Y are  available  for  every  5 days  covering  a | 

month.  These  values  are  both  rj^  and  smoothed.  They  are  available 

one  month  late.  These  values  are  obtained  from  both  satellite  [ 

tracking  and  astronomic  latitude  data.  i 

3.  DPMS  (Doppler  Polar  Monitoring  Service).  Values  of  X and  Y,  | 

available  every  5 days,  approximately  two  weeks  late.  These  values 

are  obtained  from  doppler  tracking  of  satellites.  Use  the  Rapid  | 

Service  time  (UTl  - UTC)  or  the  BIH  raw  time  with  these  values.  | 
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4.  IPMS  (International  Polar  Motion  Service).  Values  of  X and  Y, 
available  about  every  18  days  (every  0.05  year  based  on  the  Besselian 
year)  These  values  are  both  final  (about  1/3  year  late)  and  preliminary 
(about  1/4  year  late).  The  values  are  based  on  a multistation  solution 
using  observ/at ions  from  about  50  astronomical  latitude  observatories. 

5.  ILS  (International  Latitude  Service).  Values  of  X and  Y,  available 
for  every  0.05  part  of  the  Besselian  year.  These  values  are  about  1/3 

year  late.  They  are  a weighted  average  of  the  values  from  6 specific  ; 

astronomical  latitude  observatories.  The  BIH  smooth  times  should  be  !i 

used  here  and  with  source  4 above.  ;* 


For  comparison  purposes.  Tables  13,  14,  and  15  show  the  various  values  at  4 times 
during  1976.  All  values  in  Tables  13  and  14  are  "UTl"  - "UTC". 

TABLE  13 

Preliminary  Times  for  Use  in  Polar  MoUon  Tests 


..Kl 

II 


Date 

Extrapolated  Estimated 

, seconds  seconds 

DUTl 

seconds 

mill 

DUTOOT* 

' da^ 

Ji 

1 

2/1 

+.626  +.636 

+.6 

-2.5 

''  ii 

■ 

4/1 

+.470  +.454 

+ .4 

-3.3 

f 

F, 

L 

7/1 

+.182  +.179 

+ .2 

-2.7 

V 

i 

10/1 

-.044  -.047 

-.1 

-3.0 

i 

i 

* 

DUTOOT  is  the  rate  of  change  of  ’UTl  - UTC" 

b 

h 

.! 

TABLE  14 

i 

Final  Times  for  Use  in  Polar  Motion  Tests 

1 

Rapid  Service  BIH  - Raw 

BIH  - 

Smooth 

3 

1 

.1 

Date 

seconds 

seconds 

seconds 

2/1 

+.6417 

+.6399 

+.6393 

4/1 

+.4587 

+.4575 

+.4570 

7/1 

+.1862 

+.1880 

+.1865 

10/1 

-.0510 

-.0518 

-.0509 
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B.  Effect  on  Look  Angles 

The  follow^.'y  comparison  runs  were  made. 


Time  Only  Time,  X,  Y 


1. 

Extrapolated 

8. 

Rapid  Service  Values 

2. 

Estimated 

9. 

BIH  Smoothed  Values 

3. 

DUTl 

10. 

BIH  Raw  Values 

4. 

Rapid  Service  Time 

11. 

DPMS,  BIH  (Raw)  Time  Values 

5. 

BIH  (Raw)  Time 

12. 

DPMS,  BIH  (Smoothed)  Time  Values 

6. 

BIH  (Smoothed)  Time 

13. 

ILS,  BIH  (Smoothed)  Time  Values 

7. 

No  Correction 

14. 

IPMS,  BIH  (Smoothed)  Time  Values 

The  first  and  most  important  comparison  will  be  between  the  uncorrected  or 
"std"  run  and  both  Rapid  Service  runs  (with  and  without  X,  Y,  values). 

The  DUTl  corrected  run  is  also  compared.  This  comparison  is  shown  in  Table  16 
in  the  sense  "Rapid  Service"  - "Std."  and  "DUTl"  - "Std.". 

TABLE  16 

Comparison  of  Polar  Motion  Corrected  Data  with  Uncorrecced  Data 


RS  (Time  only) 

RS  (Time 

. X,  Y) 

DUTl 

AZ 

EL 

AZ 

EL 

AZ 

EL 

Time 

arcseconds 

arcseconds 

arcseconds 

1:00-1/1 

-.25 

-.30 

-.60 

-.30 

-.50 

1 

OY 

CD 

4:00 

0. 

-.40 

-.40 

1 

o 

-.04 

-.83 

7:00 

+ .83 

-.25 

+ .43 

-.25 

+1.73 

-.50 

10:00 

+.50 

-.30 

+ .11 

-.30 

+1.04 

+ .04 

13:00 

-.11 

+ .40 

-.50 

+ .40 

-.22 

+ .83 

16:00 

-.29 

+ .29 

1 

CD 

+ .29 

-.56 

+ .61 

19:00 

-.40 

+ .11 

-.79 

+ .11 

-.83 

+ .18 

22:00 

-.40 

-.14 

-.79 

-.14 

-.79 

-.29 

23:30 

-.30 

-.25 

-.72 

-.25 

-.68 

I 

o 

As  can  be  seen,  the  DUTl  correction  may  overcorrect  as  much  as  1.5  arcseconds. 
In  runs  such  as  these  the  Estimated  correction  is  the  best  long-term  correction 
available. 
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Let  us  compare  the  RS  (Time  Only)  run  aga.nsc  the  other  time  runs  (1,  2,  b, 
and  6).  There  will  be  little  difference  in  the  results,  tabulated  in  Table  17 
in  the  sense  “hun"  - "RS".  The  results  of  the  other  runs  with  X,  Y data  (numbers 
9 to  14)  will  also  be  compared  with  the  "RS  with  X,  Y oata"  run.  This  comparison 
is  listed  in  Table  18  in  the  same  sense  ("Run"  - "RS").  Again,  there  is  little 
difference  in  the  effect  of  the  uvi'ferent  data.  Finally,  tne  data  from  runs 
with  X,  Y data  will  be  differenced  with  data  from  runs  with  only  time  corrections 
(the  same  as  used  for  the  X,  Y run).  These  differences  are  listed  in  Table  19 
in  the  sense  "Run  with  X,  Y"  - "Run  with  time  only".  These  differences  are  in 
the  azimuth  channel  only.  Note  that  the  difference  simulates  a simple  bias. 

The  input  values  for  time  errors  for  the  various  runs  are  listed  below. 

Values  Used 
DUTDOT  = -3.7  milliseconds/day 


Type  Extrap  Estim.  DUTl 

Rpd.  Srv. 

D.  BIH 
Raw 

Smooth 

DPMS 

IPMS 

ILS 

Time  (seconds)  -.041  -.04037  -.1 

-.0476 

-.04876 

-.04766 

- 

- 

— 

X (arcseconds)  - - _ 

+ .272 

+ .217 

1-.2448 

^.1408 

+ .215 

+ .145 

Y (arcseconds)  _ - - 

+ .161 

+.1526 

1- . 1 606 

+.1538 

+ .172 

+.145 

TABLE  17 


Comparison  of 

the  Effects  of 

Time 

Estinia  tes 

of  Polar 

Motion 

Extrap  Time  Estimated 

Time 

BIH  Time  (Raw) 

BIH  Time 

(Smooth) 

DUTl 

Date 

AZ 

EL 

AZ 

EL 

AZ 

EL 

AZ 

EL 

AZ 

EL 

UTC 

arcseconds 

arcseconds 

arcs 

econds 

arcseconds 

arcseconds 

1:00-1/1 

-.04 

+ .04 

+ .04 

+ .07 

0. 

0. 

0. 

0. 

-.25 

-.36 

4:00 

0. 

+ .04 

0. 

+ .04 

0. 

0. 

G. 

0. 

-.04 

-.43 

7:00 

-.11 

+ .04 

-.14 

+.04 

0. 

0. 

0. 

0. 

+ .90 

-.25 

10:00 

-.07 

-.04 

-.07 

-.04 

0. 

0. 

0. 

0. 

+ .54 

+ .36 

13:00 

+.04 

-.07 

+.04 

-.07 

0. 

0. 

0. 

0. 

-.11 

+ .43 

16:00 

+ .04 

-.04 

+ .07 

-.04 

0. 

0. 

C. 

0. 

-.29 

+ .30 

19:00 

+.04 

-.04 

+ .07 

-.04 

0. 

0. 

0. 

0. 

-.43 

+ .07 

22:00 

+ .07 

+ .04 

+ .07 

+.04 

0. 

0. 

0. 

0. 

-.40 

-.14 

23:30 

+.04 

+ .04 

+ .04 

+.07 

0. 

0. 

0. 

0. 

-.36 

-.25 

Only  the  DUTl  correction  produced  large  differences  from  tne  Rapid  Service 
post-epoch  data.  For  precision  work  the  DUTl  correction  can  not  be  used. 
Note  also,  from  Table  16,  that  failure  to  correct  polar  motion  causes 
relatively  large  errors. 
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TABLE  18 

Comparison  of  the  Effects  of  Time  and  Polar  Coordinate  Estimates  of  Polar  Motion 


C.  Conments 


In  view  of  Table  16,  DUTl  should  not  be  used  In  this  program.  For  best 
results  use  the  Extrapolated  values  for  the  day  in  question.  The  Estimated 
value  (good  for  60  days  in  advance)  will  also  produce  good  results. 

For  post-epoch  work  the  Rapid  Service  values  are  as  good  as  any.  However, 
in  view  of  the  differences  shown  in  Table  14  for  different  epochs,  the 
Rapid  Service  values  must  be  checked  for  the  occurrence  of  large  differences 
between  these  values  and  the  refined  BIH  raw  values.  The  smoothed  BIH  values 
might  be  better  for  some  applications  such  as  long-term  studies. 

The  various  values  of  X,  Y available,  with  the  exception  of  the  ILS  values 
(which  should  be  used  only  with  studies  requiring  data  available  before  1965), 
are  almost  equal.  In  particular,  the  difference  between  any  two  modulii 

V2  2 

(x  + y ))  is  on  the  order  of  .05  arcseconds.  Due  to  the  interpolation 
required,  the  IPMS  values  are  less  accurate  than  the  BIH,  or  DPMS,  or  Rapid 
Service  data  near  the  middle  of  an  interpolation  period. 

4.7  Refraction  Correction 

Provision  has  been  made  in  STARl  to  accommodate  different  refraction  subroutines. 
However,  only  a simple  algorithm  (Refrac  Method  #1)  has  been  implemented.  In 
Table  20,  a typical  example  is  shown.  As  usual.  Delta  = EL  (refracted)  - 
EL(true).  EL(true)  is  also  listed.  There  is  no  azimuth  correction  in  this 
model.  The  refraction  equation  is  given  below,  followed  by  Table  20. 

EL(refracted)  = EL(true)  + 2.73D-4*RAD/DTAN(EL(true)) 

The  above  equation  is  good  for  optical  refraction.  For  RF  refraction, 
replace  "2.73"  with  "3.36". 


^ The  experience  at  WTR  with  this  equation  indicates  it  is  surprisingly 
accurate  at  elevations  above  15**.  The  two  numerical  values  can  be  computed 
from  formulas  (for  N x 10“®)  on  pages  2 and  3 of:  Landry,  P.M.  and  Parks, 
L.D.;  "Atmospheric  Refraction  Effects  on  Tracking  System  Data";  APGC  Tech. 
Rpt.  #APGC-TDR-63-28  (1963).  Input  values  are:  T * 62®F,  e » 15  mb, 

P = 1008  mb  (effective  height  of  250  meters). 
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TABLE  20 

Comparison  of  Refraction  Corrected  Data  with  Uncorrected  Data 


EL  (True) 

Delta 

Time 

degrees 

arcminutes .arcseconds 

1:00-1/1 

41.11157 

+1'  4.53" 

2:00 

33.04827 

1'  25.09" 

3:00 

25.50169 

1*  58.06" 

4:00 

18.72165 

2'  46.14" 

5:00 

12.95303 

4'  4.82" 

6:00 

8.43667 

6'  19.63" 

7:00 

5.39193 

9'  56.61" 

8:00 

3.98658 

13'  28.00" 

9:00 

4.30438 

12'  28.13" 

10:00 

6.32602 

8'  27.93" 

11:00 

9.93365 

5'  21.55" 

12:00 

14.93664 

3'  31.07" 

13:00 

21.10474 

2'  25.87" 

14:00 

28.19450 

r 45.04" 

15:00 

35.96043 

r 17.60" 

16:00 

44.14615 

58.03" 

17:00 

52.44546 

43.32" 

18:00 

60.39787 

31.37" 

19:00 

67.12341 

23.75" 

20:00 

70.87796 

19.80" 

21:00 

69.76861 

20.74" 

22:00 

64.45318 

26.91" 

23:00 

57.04342 

36.51 " 

4.8  Change  In  Apparent  Position 

The  final  result  treated  in  this  section  Is  the  change  In  apparent  star 
position  over  a short  time  period.  Table  21  below  has  a comparison  of 
apparent  star  positions  of  the  same  star  over  1/2  hour  intervals.  It  Is 
clear  from  this  table  that,  for  high  accuracy,  the  apparent  position  must 
be  updated  to  the  time  of  look  angle  computation  rather  than  dally. 
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TABLE  21 

Effect  of  Apparent  Position  Updating 


Time 

Declination 

Right  Ascension 

Degrees 

Min 

Seconds 

Hours 

Min 

Seconds 

1:00 

-8 

13 

21.9350 

+5 

13 

23.1345 

1:30 

-8 

13 

21.9330 

+5 

13 

23.1342 

2:00 

-8 

13 

21.9310 

+5 

13 

23.1338 

2:30 

-8 

13 

21.9291 

+5 

13 

23.1335 

3:00 

-8 

13 

21.9271 

+5 

13 

23.1331 

3:30 

-8 

13 

21.9252 

+5 

13 

23.1328 

Declination  (1:00)  - Declination  (3:30)  = 0.0098  arcseconds 

RT.  Ascension  (1:00)  - RT.  Ascension  (3:30)  = 0.0017  seconds  = 0.0255  arcseconds 


VERIFICATION  OF  SUBROUTINES 


5.0  INTRODUCTION 

The  subroutines  used  In  STARl  were  verified  by  direct  calculation  and  comparison 
for  the  simple  subroutines  MV,  TIMEE,  OCT,  OCTO,  and  ASDCRK.  Simple  debug 
runs  of  both  TIMEE  and  MV  are  listed  below.  ASDCRK,  OCT,  and  OCTO  require 
special  consideration.  The  remaining  routines  were  previously  verified 
(VERNAL),  required  no  verification  (BLOCK  DATA),  or  (for  CNSTNT  or  CNSTJ) 
will  be  verified  by  comparison  with  data  from  reference  6 (hereafter  called 
AENA) . 


5.1  Verification  of  MV  and  TIMEE 


MV  multiples  the  matrices  Q and  V to  obtain  the  matrix  0.  The  matrix  Q is 
a 3x3  matrix  and  V and  0 are  3x1  vectors.  In  this  case. 


V 


0.19995884687 

0.96921698587 

-0.14315545789 


and 


Q = 


1 . 

1.1271198623  x IC'^ 
4.8937153183  x lO'^ 


-1.1271198623  x 10  ^ 
1 . 

-2.9917810259  x 10'^ 


-4.8937153183  x 10 
2.9917810259  x 10 
1 . 


The  results  are; 


Calc. 

Diff. 

0(1) 

= 0.1998566106 

0.1998566101 

-5  X 10‘^° 

0(2) 

= 0.9692352415 

0.9692352408 

o 

1 

o 

X 

t 

0(3) 

= -0.143176693 

-0.1431746693 

0. 

where  DIFF  = Calc.  - MV,  The  difference  is  due  to  roundoff  and  lack  of 
precision  in  the  calculated  portion.  The  MV  column  above  is  taken  from 
the  computer  printout. 


57 


We  win  similarly  verify  TIMEE.  The  TIMEE  column  is  taken  from  a computer 
printout.  Below,  the  difference  between  GMT  and  local  time  is  +8  hours. 


r' 


i 


Variable 

TIMEE 

Calc. 

LD 

2 

2 

LS 

3600 

3600 

IDG 

2 

2 

GH 

1.0 

3600/3600  ^ 

IHG 

1 

1 

MG 

0. 

1.0  - 1 = b 

IDL 

2 

2 

XLH 

-7 

1 

00 

H 

1 

IDL 

1 

2 - 1 » 1 

XLH 

17.0 

-7.  + 24.  = 

IHL 

17 

17 

ML 

0. 

17.0  - 17  =>  0. 

Hence,  TIMEE  and  MV  are  verified. 

5.2  Verification  of  ASDCRK 

The  ASDCRK  subroutine  contains  both  of  the  transformati Ohi  X,  Y,  Z cosines 
to  Right  Ascension  - Declination  coordinates;  that  is,  (X,  Y,  Z)  to  (RAS, 

DEC)  and  the  converse  (RAS,  DEC)  to  (X,  Y,  Z).  Hence,  the  double  trans- 
formation (RAS,  DEC)  to  (X,  Y,  Z)  and  (X,  Y,  7)  to  (RAS,  DEC)  was  repeated 
20  times  for  example  verification.  The  final  values  of  (RAS,  DEC)  were 
then  compared  with  the  initial  values.  The  results  were; 

RAS(f)  - RAS(i)  < 5.  X 10‘^^  radians 

DEC(f)  - DEC(i)  = -1.  X 10‘^®  radians  = -2.05  x lO"^^  arcseconds 

The  initial  and  final  X,  Y,  Z components  could  also  be  compared.  These  results 
were: 

X(f)  - X(i)  < 5.  X 10'^^  radians 

Y(f)  - Y(i)  < 5.  X 10“^^  radians 

Z(f)  - Z(i)  < 5.  X 10"’^  radians 
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Further, 


|X^  + - l.|  <2.  X 10’^^ 

One  can  conclude  that  one  part  of  the  transformation  js,  the  inverse  of  the 
other  part.  So,  all  that  one  must  do  is  verify  one  side  of  this  subroutine. 
The  side  chosen  is  the  transformation  from  (X,  Y,  Z)  to  (RAS,  DEC),  since 
this  involves  the  arctrig  functions. 

Initially,  the  transformation  was  written: 


i 


DEC  =>  DARSIN(X(3)) 


where 


X(l)  = X.  X(2)  = Y,  and  X(3)  = Z. 

This  result  was  compared  with  the  1^-  place  trig  subroutines  on  the  WANG. 

The  results  differed  by  as  much  as  20  arcseconds. 

By  referring  to  the  IBM  manual  of  system  subroutines,  it  was  discovered  that 
the  DATAN  subroutine  was  more  accurate  in  the  small  angle  region.  Jerry  Trimble 
of  FEC/ Performance  Analysis  Department  provided  the  author  unpublished 
results  of  a numerical  comparison  of  the  IBM  DATAN  and  DARSIN  subroutines. 

DATAN  was  as  much  as  10  times  more  accurate  for  small  angles.  Hence,  the 
transformation  was  changed  to: 


ii 

’j 


DEL  = DSQRT(X(2)*X(2)  + X(1)*X(1)) 
DEC  = DATAN(X(3)/DEL) 

The  other  coordinate  transformation  is: 


HAS  = DATAN2(X(2),X(1)) 

These  transformations  were  checked  with  many  values  of  (X,  Y,  Z)  to  prove 
that  the  resulting  angles  were  in  the  correct  quadrant.  The  numerical  values 


b9 


1 


1 

i 

■i 

,1 


i 


were  compared  with  the  WANG  1^-place  trig  pack  values.  The  results  are 
listed  below: 

DEL(STARl)  - DEL(WANG)  < 4 x 10'^^  = 8 x 10'^  arcseconds , 

DEC(STARl)  - DEC(WANG)  < 3 x 10"^^  radians  ■ 6 x 10”^  arcseconds,  and 

RAS(STARl)  - RAS(WANG)  < 9 x 10"^^  radians  * 2 x 10'^  arcseconds. 

These  results  are  at  the  limit  of  the  WANG's  accuracy. 

Selected  values  were  then  computed  by  hand  using  the  tables  available  In 
reference  11  (see  page  92  and  Table  4.14).  The  maximum  difference  was  less 
than  0.0005  arcseconds  between  values  computed  from  the  tables  and  values 
computed  via  ASOCRK.  Still,  this  subroutine  will  fail  for  declinations  of 
90°  and  -90°.  Recall  that  this  program  presupposes  that  no  star  has  such 
declination.  The  Right  Ascension  Calculation  has  no  such  problem. 

5.3  Verification  of  OCT  and  OCTO 

These  two  subroutines  can  be  verified  simultaneously  by  comparing  octal  output 
against  angle  Input  for  several  encoder  bit  values.  This  Is  true  because  OCTO 
sets  the  parameters  for  OCT.  The  only  part  of  OCT  not  verified  in  this  manner 
is  the  round-up  loop.  Examination  of  the  example  In  Section  3 will  provide 
verification  of  this  part  of  the  program  except  for  the  values  of  LIP  provided 
by  OCTO.  A simple  hand  calculation  was  sufficient  to  verify  these  three  values. 

The  verification  of  the  central  portion  of  OCT  (and  OCTO)  proceded  in  two 
steps.  First,  the  octal  number  was  changed  to  base  10;  second,  this  number 
was  multiplied  by  the  LSB  (least  significant  bit  size)  of  that  encoder.  The 
resulting  angle  was  then  compared  with  the  Input  angle  - specifically,  DELTA 
= Input  Angle  - Calculated  Angle.  These  results  are  summarized  in  Table  23. 

A list  of  LSB  versus  encoder  bit  size  can  be  found  In  Table  22. 
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TABLE  22 


Part  1 


TABLE  OF  ANGLE  LSB  VALUES 


DEGREES  BIT 

180  1 

90  2 

45  3 

22.5  4 

11.25  5 

5.625  6 

2.8125  7 

1.40625  8 

0.70312b  9 

0.3515625  10 

0.17578125  11 

0.087890625  12 

0.0439453125  13 

0.02197265025  14 

0.010986328125  15 

0.0054931640625  16 

0.00274658203125  17 

0.001373291015625  18 

0.0006866455078125  19 

0.00034332275390625  20 

0.000171661376953125  21 

0.0000858306884765625  22 

0.00004291534423828125  23 

0.000021457672119140625  24 

0.0000107288360595703125  25 

0.00000536441802978615625  26 

0.000002682209014892578125  27 

0.0000013411045074462890625  28 

0.00000067055225372314453125  29 

0.000000335276126861572265625  30 

0.0000001676380634307861328125  31 

0.00000008381 9031 71 539306640625  32 

0 . 00000004 19095 1 5857696533203125  33 


MILS 


3200. 

1600. 

800. 

400. 

200. 

100. 

50. 

25. 

12. 

5 

6. 

25 

3. 

125 

1. 

5625 

.78125 

. 390625 

.1953125 

.09765625 

.048828125 

.0244140625 

.01220703125 

.006103515625 

.0030517578125 

.00152587890625 

.000762939453125 

.0003814697265625 

.00019073486328125 

.000095367431640625 

.0000476837158203125 

.00002384185791015625 

.000011920928955078125 

.0000059604644775390625 

.00000298023223876953125 

.000001490116119384765625 

. 00000074505805969238281 25 

TOTAL 


359  99999995 '090484142303466796875 
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6399.9999992549419403076171875 


Kart  z 


ARCSECONDS 

BIT 

MILLIRADIANS 

648000. 

1 

31 41 . 592653589793238462644 

324000. 

2 

1 570. 79632679489661 9231 322 

162000. 

3 

785 . 3981 6339744830961 5661 

81000. 

4 

392 . 699081 698724154807830 

40500. 

5 

1 96 . 34954084936207740391 5 

20250. 

6 

98. 1 74770424681 038701 958 

10125. 

7 

49.087385212340519350979 

5062.5 

8 

24 . 5436926061 70259675490 

2531.25 

9 

12.271846303085129837745 

1265.625 

10 

6.135923151542564918872 

632.8125 

11 

3.067961575771282459436 

316.40625 

12 

1.533980787885641229718 

158.203125 

13 

.766990393942820614859 

79.1015625 

14 

.383495196971410307430 

39.55078125 

15 

.191747598485705153715 

19.775390625 

16 

.095873799242852576858 

9.8876953125 

17 

.047936899621426288429 

4.94384765625 

18 

.023968449810713144214 

2.471923828125 

19 

.011984224905356572107 

1.2359619140625 

20 

.005992112452678286054 

.61798095703125 

21 

.002996056226339143027 

.308990478515625 

22 

.001498028113169571514 

.1544952392578125 

23 

.000749014056584785757 

.07724761962890625 

24 

.000374507028292392878 

.038623809814453125 

25 

.000187253514146196439 

.0193129049072265625 

26 

. 000093626757073098220 

.00965645245361328125 

27 

.000046813378536549110 

.004828226226806640625 

28 

.000023406689268274555 

.002414113113403320312 

29 

.000011703344634137278 

.001207056556701660156 

30 

.000005851672317068640 

.000603528278350830078 

31 

.000002925836158534320 

.000301764139175415039 

32 

.000001462918079267160 

.000150882069587707520 

33 

.000000731459039633580 

1 295999 . 999851 1 1 793041 2292980 

% 

6283.185306448127437291754 
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B1 

Si 

IT 

ize 

1 

LSB 

deg 

1 

13 

4.39453125  x 10"^ 

17 

2.74658203  x 10'^ 

27 

2.68220901  x lO'^ 

13 

4.39453125  x 10'^ 

17 

2.74658203  x 10'^ 

27 

2.68220901  x 10'^ 

TABLE  23 

Comparisons  of  OCT  Output  with  Input 


INPUT 

deg 


J3.01368  16632  7578 

Jl. 03715  353317  120527 

n. 03729  i 726636000  1234196 

5.57545  177  127 

5.68306  4025  2069 

1.53020  2132203  570499 


-14.03 

-.09 

-0.51 

-.05 

-.007 

-.68 

-20.18 

-.13 

+1.37 

+ .14 

+ .009 
1 

+ .91 

ARS  = arcsecond 

tt  Isb  = fraction  of  the  LSB  value 

All  results  are  within  an  LSB  of  the  input  with  respect  to  the  particular 
encoder  bit  size. 


5.4  Verification  of  CNSTNT  and  CNSTJ 

A.  CNSTJ  ENTRY 

The  subroutine  CNSTJ  was  verified  by  using  it  to  compute  its  output  values,  J 
and  J' , at  ten  day  intervals  throughout  1976  and  comparing  the  output  against 
AENA  (pages  328  to  331).  In  1924  entries  there  were  only  28  errors  - of  a 
maximum  value  of  1 unit  (In  J the  units  are  1 x lO”^  seconds;  in  J'  they  are 
1 X 10’^  arcseconds).  The  errors  were  equally  divided  between  the  northern  and 
soutnern  hemispheres.  The  average  error  (sum  of  errors  divided  by  28)  is 
-1.8  X lO"^  arcseconds  or  approximately  -0.2  units.  Hence,  the  errors  are  thought 
to  be  due  to  rounding  only. 
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B. 


Constants  Computed  in  CNSTNT 


The  subroutine  CNSTNT  required  more  involved  verification.  The  constants 
computed  in  CNSTNT  are  E (Earth's  obliquity),  EOFE  and  DEFE  (difference 
between  the  "Apparent"  and  "Mean"  Vernal  Equinox  and  the  first  derivative 
of  that  difference),  and  DW  (secular  change  in  Earth's  spin  rate).  Also, 
TTAU  (fraction  of  a tropical  year  since  the  nearest  Bessel ian  New  Year)  is 
passed  via  Comnon.  Finally,  one  must  also  compute  DB  (the  UTC  date  of  the 
Bessel ian  New  Year),  DL  (the  length  of  that  year),  and  XM  (the  annual  pre- 
cession in  declination). 

First,  TTAU  was  compared  against  AENA  (pages  308  to  322)  for  a full  year. 

No  difference  was  found  for  four  decimal  places  (the  accuracy  of  the  numbers 
listed  in  AENA).  UL  and  DW  are  single  line  computations  of  small  errors. 
They  were  verified  by  hand  calculation.  The  equations  can  be  found  in 
reference  2 (pages  JIO  and  J8).  Typical  examples  are: 

-4 

DW  * -3.26  X 10  radians  per  second,  and 
DL  = 365.242198698  days,  both  for  year  1976. 

The  constant  DB  has  already  been  partially  verified,  as  it  must  be  known 
in  order  to  compute  TTAU.  However,  a comparison  table  against  values  found 
in  reference  4 (pages  434  - 435)  is  listed  below. 


TABLE  24 

Bessel ian  New  Year  Data  Test 


Year 

Calculated  Date 

AENA  Date 

1900 

Jan  0.81352 

Jan  0.813 

1964 

Jan  1.3142075 

Jan  1.314 

1975 

Jan  0.9783811 

Jan  0.978 

1976 

Jan  1.2205786 

Jan  1.221 

1977 

Jan  0.4627761 

Jan  0.463 
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Two  values  of  E are  compared  with  AENA  values  below.  The  difference  Is  due 
to  lack  of  precision  In  the  calculator  used  to  change  radians  to  degrees, 
minutes,  and  seconds. 


E = 23“26'32.656"  (AENA)  or  23‘’26'32.655"  (calc.)  for  year  1976 
E = 23°26'32.187"  (AENA)  or  23°26'32.186"  (calc.)  for  year  1977 

Note  that  the  equation  for  EOFE  is:  (AA(1)  - XM)/ETN,  where  XM  is  the  annual 
precession  in  declination.  Also,  XM  is  approximately  20  arcseconds/year 
and  is  the  same  numerical  size  as  AA(1).  Hence,  if  the  computation  for  EOFE 
can  be  verified,  the  calculation  of  XM  will  also  be  verified.  Since  EOFE  + DEFE  = 
EOFE  for  the  next  day,  one  can  easily  check  DEFE  also. 

Below  is  a comparison  of  calculated  values  of  EOFE  with  those  available  in 
AENA  (pages  12  - 19). 


TABLE  25 

Equation  of  Equinox  Comparison,  Year  1976 


Note  that  the  values  agree  to  the  nearest  decimal  that  AENA  lists.  This 
completes  the  discussion  of  the  constants  computed. 
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C.  First  and  Second  Derivatives  of  Bessel ian  Numbers 

The  assignment  of  Block  Data  values  to  the  proper  day,  etc.,  was  checked  by 
hand  against  the  printout.  The  derivatives  are  computed  by  using,  for  the 
first  derivative,  a simple  forward  difference  (for  DAP,  a backward  difference). 

For  the  second  derivative  a symmetric  second  difference  is  used.  These  techniques 
were  checked  by  hand.  They  are  valid  if  the  differences  are  smaller  than  the 
approximated  function  in  decreasing  order  and  if  they  change  slowly.  A comparison 
of  many  values  of  the  function  and  first  and  second  derivatives  revealed  that 
the  value  of  the  first  derivative  is  usually  1/10  or  less  of  the  value  of  the 
initial  function,  while  the  value  of  the  second  derivative  averages  1/2  to  1/10 
the  value  of  the  first  derivative. 

Similar  graphs  of  the  changing  values  of  the  numbers  A,  B,  C,  D,  and  f, 
together  with  their  first  and  second  derivatives,  are  shown  in  Figures  4,  5, 

6,  7,  and  8.  For  all  of  the  graphs  "1"  denotes  a graph  of  the  Besselian  Day 
Number,  "2"  denotes  a graph  of  the  Day  Number's  1st  differences,  and  "3"  denotes 
the  graph  of  the  2nd  differences.  Note  that  only  the  absolute  value  of  the 
parameters  was  graphed.  The  graphs  were  prepared  for  the  first  30  days  of 
1976,  but  are  representative  of  changes  throughout  the  year.  The  units  are 
arcseconds  for  A,  B,  C,  and  D,  but  are  seconds  of  time  for  f.  Except  for 
Figure  5 (the  number  B),  the  previous  discussion  holds.  For  B the  first  and 
second  derivatives  can  be  equal,  and  the  second  derivative  can  be  larger  than 
the  first  derivative.  However,  both  derivatives  are  100  times  smaller  than 
8 itself. 

5.5  Verification  of  VERNAL 

VERNAL  is  a subroutine  taken  from  reference  8^  (page  26).  The  verification  of 
this  subroutine  can  be  found  in  the  same  reference  on  pages  34  and  37.  In  order 
to  improve  this  subroutine  even  further,  the  initial  (1973)  value  was  recomputed 
by  using  the  equations  of  reference  2 (page  42).  Examination  of  the  second  order 
terms  in  references  2 and  8 showed  that  the  time  in  the  second  order  term  required 
a small  change.  Also,  the  leap  year  update  was  changed  from  mixed  mode  arithmetic 
to  integer  arithmetic  in  order  to  avoid  ambiguity.  These  changes  are  listed  below. 

t The  factor  of  -1/2  on  page  26  is  a misprint  and  should  be  eliminated. 
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FIGURE  • INDEPENDENT  DAY  NUMBER  f AND 

1ST.  AND  2ND.  DERIVATIVES 


Examination  of  differences  between  the  two  forms  of  VERNAL  for  year  1976 
showed  that  change  1.  had  no  observable  effect,  change  3.  had  an  effect  of 
about  0.095  arcseconds  at  the  equator,  and  change  2.  has  an  effect  of  about 
0.005  arcseconds.  For  VAFB  stations  the  maximum  total  effect  was  on  the  order 
of  0.04  arcseconds. 

5.6  Discussion  of  BLOCK  DATA 

The  primary  data  contained  In  BLOCK  DATA  are  the  Bessel Ian  Day  Numbers  A,  B,  C, 
and  D and  the  Independent  Day  Number  f.  These  were  described  In  Section  2. 

The  other  data  contained  In  BLOCK  DATA  are  Initial  values  for  other  parameters 
and  basic  constants.  This  data  can  be  checked  by  Inspection. 


However,  the  Bessel  Ian  and  Independent  Day  Numbers  must  be  verified  by  some 
other  means.  The  test  program  used  to  produce  the  printout  for  the  verifi- 
cation of  CNSTJ  provides  such  a means.  All  of  the  Bessel  Ian  Day  Numbers 
are  used  to  calculate  J and  J'.  Hence,  a run  printing  out  every  day  of  the 
year  can  be  used  to  verify  A,  B,  C,  and  D by  comparing  against  the  AENA 
tables.  The  values  for  days  In  between  the  10-day  values  In  AENA  can  be 
checked  by  Interpolation.  Refer  to  Section  3 for  the  equations  used. 

It  is  known  that: 

f ■ m T + EOFE, 

where 

m - 3.07234  + 1 .86D-3*T/36524.22, 

T ■ number  of  days  from  1900.0  to  the  nearest  Bessel  Ian  New  Year, 
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r 

T = fraction  of  a tropical  year  from  date  to  the  nearest  Bessel  Ian 
New  Year,  and 

EOFE  = Equation  of  Equinoxes. 

Hence,  EOFE  = f - nrc . 

A computer  program  called  JJJPP  has  been  designed  to  check  the  Bessel ian  Day 
Numbers  and  f for  a given  year  in  the  manner  described  above.  This  was  accom- 
plished for  the  BLOCK  DATA  for  1976.  Each  year  the  new  BLOCK  DATA  can  be  so 
verified.  The  program  JJJPP  is  listed  in  Appendix  VIII.  This  completes  the 
verification  of  the  subroutines. 
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6.0 


VERIFICATION  OF  MAIN 


This  program  was  verified  by  comparison  against  certified  WTR  programs. 

Of  course,  this  program  should  be  more  accurate  than  previous  programs  at 
SAMTEC,  so  verification  does  not  mean  absolute  agreement.  The  programs 
used  for  verification  were  the  "mean"  to  "apparent"  transformations  on  the 
Sigma  5 computer  at  the  FPQ-14  and  the  "look  angle"  output  of  STASHO.  The 
version  of  STASHO  used  was  that  available  on  the  WANG  at  PAD.  Verbal  communi- 
cation from  Bob  Baker  of  FEC/ Performance  Analysis  Department  informed  the 
author  of  an  unpublished  comparison  between  the  WANG  version  of  STASHO  and 
the  7094  version  which  showed  agreement  between  these  programs  within  the 
precision  of  the  program's  printout  (0.2  arcseconds).  The  two  remaining 
effects  needing  verification  were  the  diurnal  aberration  effect,  which  was 
checked  against  Table  2.6  of  reference  4,  and  the  polar  motion  correction. 

This  last  correction  was  verified  by  hand  calculation. 

6.1  Verification  of  the  "Apparent"  Position  Computation 


"Apparent"  star  positions  at  0.0  hour,  UTC,  for  the  dates  listed  in  Table  P7 
were  compared  against  those  computed  by  the  FPQ-14  computer.  The  FPQ-14  data 
is  the  4th  Fundamental  Catalog  of  1950.0  (FK4).  For  position  (Right  Ascension, 
declination)  the  Mean  Star  Table  of  the  AENA  for  the  appropriate  year  (1976 
before  July  1,  1977  afterward)  was  used.  For  both  proper  motion  and  parallax 
the  data  in  reference  9 was  used.  By  use  of  radial  velocities  obtained  from 


reference  10  and  the  data  in  reference  9,  the  proper  motion  parameters  were 


Mean  Star  Coordinates 


A SCO 


Note  that  the  initial  values  are  accurate  to  only  +.5  arcsecond  in  declination 
and  +.05  second  in  Right  Ascension.  This  is  much  lower  than  the  accuracy  of 
the  4th  Fundamental  Catalogue.  Nevertheless,  in  Table  27  below  good  agreement 
is  obtained  between  the  two  programs.  The  differences  are  tabulated  in  the 
usual  manner;  "STARl"  data  - "FPQ-14"  data.  The  average  error  in  Table  23  is 
-.214  arcseconds  in  declination  and  +.018  seconds  in  Right  Ascension. 


TABLE  27 

Comparison  of  Apparent  Positions  Obtained  from 
STARl  Data  and  from  FPQ-14  Sigma  5 Data 


Star  Name 

Date 

Right  Ascension 

Decl ination 

UTC 

seconds 

arcseconds 

Arcturas 

1/5/1976 

-.006 

-.402 

Betel guese 

II 

+ .026 

-.606 

Castor 

II 

+.036 

-.313 

Dubhe 

It 

-.017 

+ .175 

A1  ioth 

5/10/1976 

+.007 

-.558 

Procyon 

II 

+.069 

-.039 

Rigel 

II 

+ .030 

-.039 

Alhena 

8/15/1976 

+ .026 

-.486 

Bellatrix 

H 

+.064 

+ .080 

Vega 

M 

+ .010 

-.212 

Aldebaran 

11/19/1976 

+ .001 

-.144 

Antares 

II 

-.006 

-.182 

Wezen 

II 

-.005 

-.052 
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6.2  Verification  of  the  Diurnal  Aberration  Correction 

In  reference  4,  P.50,  Is  a table  of  diurnal  aberration  corrections  for  use 
when  a star  Is  In  transit  - that  is,  when  the  local  hour  angle  is  zero  or 
twelve.  Comparisons  were  made  with  STARl  by  use  of  special  printouts  and 
linear  interpolation. 

The  apparent  star  declination  is  89.1569"  (for  Polaris).  The  site  latitude 
is  37.5®.  From  the  table  referenced  above,  the  correct  values  (after  linear 
Interpolation)  are: 

aRA  = + 17.2374,  aDEC  = 0.0,  both  in  arcseconds. 

The  values  obtained  from  STARl  are  listed  in  Fable  28.  !' 

:i 

TABLE  28  j 

Diurnal  Correction  Data  from  STARl  \ 


Set 

Local  Hour  Angle 

ARA 

ADEC 

# 

hours 

arcseconds 

arcseconds 

1 

1 1 . 7280 

-17.179 

+.0180 

1 

12.2294 

-17.191 

-.0150 

2 

23.6136 

+17.166 

-.0256 

2 

0.1150 

+17.246 

+.0071 

3 

23.5480 

+17.131 

-.0299 

3 

0.0494 

+17.251 

+.0033 
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In  Table  29  the  values  at  local  hour  angle  12.  or  0.  (obtained  by  linear 
intepolatlon)  are  listed,  together  with  their  differences  from  the  tabular 
values  (listed  In  the  form  "STARl"  - "Table").  The  results  are  within  the 
error  resulting  from  linear  Interpolation. 


TABLE  29 

Diurnal  Correction  Comparison  Test 


Set 

LHA 

aRA 

ADEC 

ARA-Table 

ADEC-Table 

# 

hours 

arcseconds 

arcseconds 

arcseconds 

arcseconds 

1 

12.0 

-17.1855 

+.0001 

-.0520 

+.00010 

2 

0.0 

+17.228 

+.00003 

-.0094 

+.00003 

3 

0.0 

+17.239 

-.00001 

+.0016 

-.00001 
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6.3 


Verification  of  "Look  Angle"  Computation 


The  data  from  STARl  was  compared  against  the  WANG  version  of  STASHO  as 
previously  mentioned.  The  Input  data  was: 

Mean  Right  Ascension  = 6 hrs. , 29  min.,  8.1  sec. 

Mean  Declination  = 56®,  24' , 22.1" 

Epoch  = 1976.0 

Site  Latitude  = +37.4958® 

Site  Longitude  = +237.4961°  E.  Long. 

All  other  parameters  (except  tne  site  encoder  size)  are  zero.  For  one  time 
period  a declination  of  36®24'22.1"  was  used.  For  all  dates  after  July  1 
the  mean  parameters  above  were  updated  to  the  1977.0  epoch  via  STAR2. 

The  comparison  results  are  listed  in  Table  30.  The  units  are  mils,  which 
were  the  units  used  by  STASHO.  Note  that  .001  mils  “ .2  arcseconds.  For 
comparison  purposes  two  STARl  decks  were  used.  The  'Full'  deck  is  the 
standard  STARl  deck.  The  'Deleted'  deck  had  all  processing  not  common  to 
both  STARl  and  STASHO  zeroed  out.  The  diurnal  aberration  correction  is  the 
largest  of  the  missing  corrections.  The  differences  in  Table  30  are  listed 
in  the  form;  "STARl"  - "WANG". 
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TABLE  30 

Comparison  of  Look  Angle  Output  Between  STAR!  and  STASHO 


Date  (UTC) 

"Deleted" 

Version 

"Full" 

Version 

Month/day 

Azimuth  Diff. 

Elevation  Diff. 

Azimuth  Diff 

. Elevation  Diff. 

Hour 

mils 

mils 

mils 

mils 

1-1/3 

+.001 

+.001 

0.0 

0.0 

7 

0.0 

+ .002 

-.001 

+.002 

13 

-.001 

-.001 

-.001 

-.001 

19 

+.001 

0.0 

0.0 

+.001 

23.5 

+.002 

+.001 

-.001 

+.001 

1-4/8 

-.003 

+ .001 

O 

o 

1 

+.002 

7 

+.001 

-.001 

0.0 

CNi 

o 

o 

1 

13 

+.001 

0.0 

+.002 

+.001 

19 

+.001 

+.002 

0.0 

+.002 
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Note  that  the  maximum  error  in  the  "Deleted"  Version  column  Is  0.6  arcseconds 
with  an  average  of  0.12  arcseconds.  Since  these  computations  are  close  to  the 
accuracy  of  the  WANG  (considering  all  computations).  STARl  is  verified. 

The  "Full"  version  has  a maximum  error  1.2  arcseconds,  which  Is  consonant 
with  the  lack  of  correction  of  the  diurnal  aberration  In  the  WANG  program. 

6.4  Termination  of  Program 

In  Section  2 the  three  abnormal  terminations  available  in  STARl  were  discussed. 
These  terminations  are  labeled  STOP'S.  They  occur  when: 

1.  The  first  day  of  the  run  is  less  than  or  equal  to  0. 

2.  The  first  day  of  the  run  preceeds  July  1 while  the  last 

day  equals  or  follows  July  1. 

3.  The  last  day  of  the  run  is  equal  to  or  greater  than 
January  2 of  next  year. 

Note  that  the  "days"  used  above  are  the  GMT  equivalents  of  the  local  start 
and  stop  days,  rather  than  the  local  day  numbers  themselves. 

Since  the  abnormal  STOP'S  are  labeled,  the  HASP  printout  will  show  this. 

Since  the  HASP  printout  is  the  first  page  of  the  run,  such  conditions  can 
be  quickly  recognized.  Both  the  HASP  and  STARl  printouts  for  the  abnormal 
cases  and  the  STARl  printout  for  the  normal  case  are  shown  in  Figure  9. 

It  is  of  interest  that  the  year  number  in  the  normal  printout  is  preset  in 

BLOCK  DATA.  Hence,  it  does  not  require  a special  change  each  year. 

6.5  Conment  on  Refraction  Corrections 

The  output  of  STARl  without  refraction  correction  is  more  than  adequate 
for  differencing  with  real  data  in  order  to  regress  for  systematic  errors 
after  the  real  data  has  been  precorrected  for  refraction  errors.  If  the 
STARl  output  'is  CO  be  used  for  pointing  instruments,  then  the  built-in 
refraction  correction  is  the  minimal  correction  needed.  However,  this 
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FIGURE  9 

Samples  of  the  Different  STARl  Terminations 
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program  was  designed  to  accommodate  more  extensive  and  elaborate  refraction 
subroutines  as  required  by  the  user.  These  subroutines  are  already  available 
at  WTR  and  need  not  be  rediscovered.  Local  experts  on  such  routines  are 
Mr.  Jerry  Trimble  and  Dr.  Ruey  Han,  both  of  FEC/Performance  Analysis  Department. 

6.6  Conclusions  and  Recommendations 

There  Is  no  doubt  that  this  program  will  produce  data  to  an  accuracy  of  0.1 
arcseconds.  Whether  or  not  the  data  is  accurate  to  0.01  arcseconds,  as  was 
the  design  criterion,  will  require  checkout  against  a more  accurate  program 
than  is  currently  available  at  SAMTEC.  The  references  used  as  background  for 
this  program  claim  an  accuracy  for  the  apparent  star  position  in  excess  of 
0.01  arcseconds. 

The  only  remaining  problem  is  that  of  changing  the  BLOCK  DATA  each  year 
(see  Appendix  V).  There  are  at  least  two  approaches  available  to  further 
simplify  this  task.  The  first  approach  would  be  to  Implement  subroutines 
to  obtain  the  additional  astronomic  data  now  missing.  If  that  happens, 
then  the  program  STARS  would  be  used  to  calculate  the  Bessel ian  Day  Numbers 
C and  D,  in  addition  to  A,  B,  and  f.  Unfortunately,  the  technical  character 
of  these  computations  requires  no  small  effort  and  time  to  program  successfully. 
The  second  approach  involves  obtaining  a tape  with  the  required  numbers  from 
the  Naval  Observatory.  Then  a program  would  have  to  be  written  to  read  the 
tape  and  output  the  data  in  useable  form. 
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USER  INSTRUCTIONS 
STARl  PROGRAM 

I 


I. 


PURPOSE 


T 


The  purpose  of  this  program  Is  to  compute  look  angles  (azimuth  and  elevation) 
for  tracking  a given  star  from  a given  site.  The  output  (for  a given  site 
and  a given  star)  is  azimuth  and  elevation  in  the  common  radar  system  - Earth 
rotating,  right  handed,  referenced  to  the  astronomic  vertical  - and  is  listed 
in  degrees  (and  decimal  degrees),  mils,  and  octal.  Other  outputs  are  time 
(both  local  and  UTC),  plunge  azimuth,  and  plunge  elevation.  The  last  two  are 
listed  in  mils  and  octal  only. 


:i 


The  octal  output  requires  special  consideration.  This  output  must  be  refer- 
enced to  the  number  of  bits  in  the  site  encoders.  Hence,  if  octal  output  is 
desired,  the  bit  size  of  the  site  encoders  must  be  known.  This  output  is 
eliminated  when  a bit  size  of  zero  is  input. 


II  Logical  Setup  of  Input  Deck 


The  following  is  a block  diagram  of  the  input  deck.  Note  that  the  first  and  last 
cards  appear  only  once  in  the  deck.  Card  blocks  2,  3,  and  5 must  appear  in  every 
run.  Card  block  4 will  appear  in  every  run  for  which  polar  motion  is  corrected, 
but  not  otherwise. 

1 Number  of  runs 


Run  Cards 


Card  for  number  of  sites  (I)  and  number  of  stars  per  site  J(I) 
for  run  1 . 

Integration  card  for  run  1 


Station  Cards 


I 


Run 


Station  coordinate  card  for  site  1 

Station  options  card  for  site  1 

3 Station  coordinate  card  for  site  2 
I Station  options  card  for  site  2 


Station  coordinate  card  for  site  I 
Station  options  card  for  site  I 
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